FINAL  REPORT 


A  Complex  Approach  to  UXO  Discrimination:  Combining 
Advanced  EMI  Forward  Models  and  Statistical  Signal  Processing 


SERDP  Project  MR-1572 


JANUARY  2012 


Fridon  Shubitidze 
Sky  Research,  Inc. 

Partner:  Thayer  School  of  Engineering,  Dartmouth 
College 


This  document  has  been  cleared  for  public  release 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

JAN  2012  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2012  to  00-00-2012 

4.  TITLE  AND  SUBTITLE 

A  Complex  Approach  to  UXO  Discrimination:  Combining  Advanced 

EMI  Forward  Models  and  Statistical  Signal  Processing 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Sky  Research,  Inc,445  Dead  Indian  Memorial  Road, Ashland, OR, 97520 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

The  research  described  in  this  report  was  conducted  in  fulfillment  of  Project  MM-1572,  ?A  Complex 
Approach  to  UXO  Discrimination:  Combining  Advanced  EMI  Forward  and  Statistical  Signal  Processing,? 
submitted  to  the  Strategic  Environmental  Research  and  Development  Program  (SERDP)  in  response  to  the 
Munitions  Management  Statement  of  Need  MMSON-07-04  ?Advanced  Technologies  for  Detection, 
Discrimination,  and  Remediation  of  Munitions  and  Explosives  of  Concern  (MEC):  UXO  Technology?.  The 
primary  objective  of  the  project  MM-1572  was  to  develop  and  validate  innovative,  robust  and  practical 
approaches  for  UXO  localization  and  classification  under  realistic  (noisy,  cluttered  background)  field 
conditions  by  combining  advanced  electromagnetic  induction  (EMI)  forward  and  statistical  signal 
processing  methodologies.  In  a  real  field  the  electromagnetic  signals  become  convoluted  with  noise  due  to 
the  instrument,  magnetic  soil  and  widespread  background  clutter.  Therefore,  the  rationale  for  a  statistical 
approach  is  to  use  an  advanced  statistical  approach  to  reduce  the  impact  of  these  noises  to  minimum.  This 
project  provides  mathematical  fundamentals,  physical  meanings  and  practical  realizations  of  forward, 
inverse  and  statistical  signal  processing  approaches  for  unexploded  ordnance  (UXO)  detection  and 
discrimination  at  live-UXO  sites.  Namely,  under  this  project  first  we  developed  and  implemented 
advanced,  physically  complete  forward  EMI  models  such  as,  the  normalized  surface  magnetic  source 
(charge/dipole)  model  (NSMS),  and  ortho-normalized  volume  magnetic  source  (ONVMS)  technique  for 
accurately  representing  the  EMI  responses  of  subsurface  metallic  targets,  then  we  combined  these 
advanced  models  with  EMI  data  inversion  approaches,  such  as  the  gradient  search  direct 
search-differential  evolution  and  etc,  for  advanced  EMI  sensor  data  inversion;  third  we  extended  the 
advanced  statistical  signal  processing  approaches,  i.e.  support  vector  machines,  Gaussian  mixture  models, 
for  discriminating  UXO  targets  from  non-hazardous  anomalies.  Finally,  the  combined  advanced  EMI 
forward  and  statistical  models  were  applied  to  ESTCP  live  site  UXO  data  sets.  Live  site  discrimination 
studies  showed  the  excellent  discrimination  performance  of  the  advanced  models  when  applied  to 
next-generation-sensor  data  collected  at  various  live  sites,  such  as  Camp  Sibert,  AL,  San  Luis  Obispo 
(SLO),  CA,  and  Camp  Butner,  NC  as  well  as  APG  test  sites.  The  technology  was  able  to  single  out  UXO 
ranging  in  caliber  from  25  mm  up  to  155  mm.  In  addition,  the  ONVMS  technique  was  seen  to  provide 
excellent  classification  in  both  single-  and  multiple-target  scenarios  when  combined  with  advanced 
multi-axis/transmitter/receiver  sensors  data. 


15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

Same  as 
Report  (SAR) 

18.  NUMBER 
OF  PAGES 

147 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


This  report  was  prepared  under  contract  to  the  Department  of  Defense  Strategic 
Environmental  Research  and  Development  Program  (SERDP).  The  publication  of  this 
report  does  not  indicate  endorsement  by  the  Department  of  Defense,  nor  should  the 
contents  be  construed  as  reflecting  the  official  policy  or  position  of  the  Department  of 
Defense.  Reference  herein  to  any  specific  commercial  product,  process,  or  service  by 
trade  name,  trademark,  manufacturer,  or  otherwise,  does  not  necessarily  constitute  or 
imply  its  endorsement,  recommendation,  or  favoring  by  the  Department  of  Defense. 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Table  of  Contents 

Table  of  Contents . ii 

List  of  Figures . v 

List  of  Tables . xi 

List  of  Acronyms . xii 

Abstract . 1 

Acknowledgements . 2 

1  Introduction . 3 

1 . 1  Background  and  objectives . 3 

1 .2  Report  structure . 4 

2  Forward  models . 5 

2.1  Introduction . 5 

2.2  The  single -dipole  approximation . 7 

2.3  Normalized  surface  magnetic  source  model . 9 

2.3. 1  Governing  equations . 9 

2.3.2  Theoretical  basis  of  NSMS . 1 1 

2.3.3  Formulation  for  bodies  of  revolution;  determining  NSMS  amplitudes  from 

data . 13 

2.3.4  The  dipole  model  as  a  limiting  case  of  NSMS . 16 

2.3.5  Invariance  properties  of  the  total  NSMS . 18 

2.3.6  Interpretation  of  the  total  NSMS . 19 

2.3.7  The  parameterized  NSMS . 19 

2.4  The  orthonormalized  volume  magnetic  source  model . 20 

2.4. 1  Orthonormal  Green  functions . 22 

2.4.2  ONVMS  procedure . 24 

2.5  Joint  diagonalization  for  multi-target  data  pre-processing . 25 

2.5. 1  The  multi-static  response  matrix . 25 

2.5.2  Interpretation  and  diagonalization  of  the  MSR  matrix . 26 

2.5.3  Algorithm  for  j  oint  diagonalization . 29 

Sky  Research,  Inc.  ii  January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


3  Inverse  Models . 31 

3.1  Introduction . 31 

3.2  Gradient-based  methods  of  optimization . 32 

3.2. 1  Stepwise  optimization . 34 

3.2.2  Simultaneous  optimization . 34 

3.2.3  Condensed  algorithm . 34 

3.3  Differential  evolution . 35 

3.4  The  HAP  method . 37 

3.4. 1  Estimating  the  location  and  orientation  of  buried  objects . 37 

3.4.2  A  simplified  HAP  method . 39 

3.4.3  Determining  the  HAP  amplitudes . 39 

3.4.4  The  HAP  method  with  gradient  information . 41 

4  Statistical  approaches . 42 

4. 1  Introduction . 42 

4.2  Mixed  modeling  applied  to  advanced  EMI  features . 43 

4.3  Gaussian  mixture  models . 48 

4.3. 1  Model-based  supervised  clustering . 48 

4.3.2  Unsupervised  classification  using  the  multivariate  normal  mixture  approach . 49 

4.4  Support  vector  machines  for  subsurface  object  classification . 5 1 

5  Advanced  models  applied  to  next -generation  sensors:  Modeling  and  validation . 54 

5 . 1  Introduction . 54 

5.2  MetalMapper . 54 

5.3  TEMTADS . 59 

5.3. 1  TEMTADS  modeling . 59 

5.3.2  APG  test-site  classification . 61 

5.4  BUD . 67 

5.5  MPV . 71 

6  ESTCP  live-site  classification  studies  using  advanced  models . 74 

6 . 1  Introduction . 74 

6.2  Camp  Sibert . 74 

6.2. 1  Target  location  and  characterization;  preliminary  pattern-matching 

classification . 75 

Sky  Research,  Inc.  iii  January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


6.2.2  SVM  classification . 79 

6.2.3  SVM  analysis  of  Camp  Sibert  data:  summary . 84 

6.2.4  Mixed  model  approach  applied  to  Camp  Sibert  data . 85 

6.3  Camp  San  Luis  Obispo  (TEMTADS,  MM.  BUD) . 87 

6.3. 1  The  total  NSMS  for  discrimination . 89 

6.3.2  SLO  discrimination  results . 89 

6.3.3  Comparisons  between  NSMS  and  Dipole  models . 94 

6.3.4  SLO  BUD  data  inversion  and  classification  studies . 97 

6.3.5  SLO  retrospective  analysis . 98 

6.4  Camp  Butner . 101 

6.4. 1  TEMTADS  data  discrimination  strategy  and  classification  results  using 

supervised  clustering . 102 

6.4.2  MetalMapper  data  discrimination  strategy  and  classification  results  using 

supervised  clustering . 106 

6.4.3  A  Comparison  between  ONVMS  and  Dipole  model . 109 

6.4.4  Camp  Butner  retrospective  study  using  semi-supervised  clustering . 110 

7  Conclusions . 119 

8  Publications . 122 

8.1  J ournal  Articles . 122 

8.2  Ph.D  Thesis . 122 

8.3  Conference  Papers . 122 

8.4  Presentations  and  posters . 124 

9  References . 125 


Sky  Research,  Inc.  iv  January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


List  of  Figures 

Figure  1:  A  dipole’s  location  in  a  global  coordinate  system . 8 

Figure  2:  The  NSMC  that  are  distributed  on  a  prolate  spheroidal  surface  is  implemented  for  a 

body  of  revolution.  The  prolate  spheroidal  coordinate  system  is  specified  by  (£  //,  <p) . 14 

Figure  3:  A  schematic  diagram  for  a  dipole  model . 17 

Figure  4:  A  metallic  object  under  the  transmitter.  The  target’s  EMI  response  at  the  receiver  coil 

can  be  calculated  from  the  equivalent  surface  or  volume  magnetic  dipole  moment  dm . 21 

Figure  5:  The  HAP  approach  for  a  dipole . 37 

Figure  6:  Determining  the  location  and  orientation  of  a  buried  target.  The  method  assumes  the 
object  is  a  point  dipole  and  exploits  an  analytic  relation  between  the  field  measured  at  r  and 

the  scalar  potential  at  the  same  point  to  find  the  location  r  ; .  The  potential  is  constructed 
using  a  layer  of  equivalent  magnetic  sources  placed  between  the  sensor  and  the  object;  r  is  a 
typical  location  on  the  layer . 40 

Figure  7:  The  MetalMapper  during  SLO  site  deployment  (left)  and  its  schematic  diagram  (right) . 55 

Figure  8.  The  MetalMapper  geometry.  The  observation  point  r  is  defined  with  respect  to  the 
global  Cartesian  coordinate  system  XYZO;  r'.  is  the  location  of  the  i-th  current  element  on 

(in  this  case)  the  T=  3  transmitter,  which  carries  a  current  /3  in  the  direction  /; . 56 

Figure  9.  Response  of  an  81-mm  mortar  illuminated  by  the  MM  Z-transmitter:  measured  (left), 

ONVMS  prediction  (center),  and  mismatch  between  modeled  and  actual  data  (right).  The 

mortar  is  placed  35  cm  below  the  sensor  center  and  oriented  45  degrees  nose  down.  The  data 

are  plotted  in  log10  scale . 57 

Figure  10.  Response  of  an  81-mm  mortar  illuminated  by  the  MM  K-transmittcr:  measured  (left), 

ONVMS  prediction  (center),  and  mismatch  between  modeled  and  actual  data  (right).  The 

mortar  is  placed  35  cm  below  the  sensor  center  and  oriented  45  degrees  nose  down.  The  data 

are  plotted  in  log10  scale . 58 

Figure  1 1 .  Response  of  an  8 1-mm  mortar  illuminated  by  the  MM  A-transmittcr:  measured  (left), 

ONVMS  prediction  (center),  and  mismatch  between  modeled  and  actual  data  (right).  The 

mortar  is  placed  35  cm  below  the  sensor  center  and  oriented  45  degrees  nose  down.  The  data 

are  plotted  in  log10  scale . 58 

Figure  12:  Photo  of  the  TEMTADS  in  deployment  at  Blossom  Point  Test  Site  (left)  and  a 

schematic  diagram  of  its  Tx/Rx  sensors  (right) . 59 

Figure  13.  Measured  (top  five  rows)  and  ONVMS-modeled  (bottom  five)  TEMTADS  data  for  a 
105-mm  projectile  at  the  25th  time  channel.  The  target  is  buried  at  a  depth  of  30  cm  and 
oriented  horizontally  relative  to  the  TEMATDS  system . 60 

Figure  14:  The  APG  TOI . 61 

Figure  15:  Comparison  between  the  inverted  and  actual  depth  for  all  65  APG  calibration  targets . 62 

Figure  16:  Inverted  total  NSMS  for  APG  test-stand  105  mm  projectile  and  81  mm  mortar . 63 

Figure  17:  Scatter  plot  of  inverted  /?  vs.  In  k  classification  features  for  APG  test-stand  TOI . 64 


Sky  Research,  Inc. 


v 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Figure  18:  Scatter  plot  of  inverted  y  vs  In  A'  (left)  and  fi  (right)  parameters  for  APG  test-stand 

TOI . 64 

Figure  19:  Scatter  plot  of  inverted  /?  vs.  In  A  classification  features  for  all  214  APG  blind-test 

anomalies . 65 

Figure  20:  Comparison  between  library  (green  lines)  and  inverted  (red  and  blue  lines)  blind-test 

total  NSMS  for  105-mm  projectiles,  81 -mm  munitions,  and  60-mm  mortars . 66 

Figure  21:  Comparisons  between  library  (green  lines)  and  inverted  (red  and  blue  lines)  blind-test 

total  NSMS  for  37-mm  and25-mm  mortars . 66 

Figure  22:  Schematic  diagram  of  the  BUD  system . 67 

Figure  23:  The  BUD  system  in  operation . 68 

Figure  24:  Comparisons  between  actual  and  predicted  data  for  an  M75  UXO  illuminated  by  the 

BUD  Z  transmitter.  Solid  lines  are  actual  data,  circles  stand  for  NSMS  predictions . 69 

Figure  25:  Comparisons  between  actual  and  predicted  data  for  an  M75  UXO  illuminated  by  the 

BUD  X  transmitter.  Solid  lines  are  actual  data,  circles  stand  for  NSMS  predictions . 69 

Figure  26:  Comparisons  between  actual  and  predicted  data  for  an  M75  UXO  illuminated  by  the 

BUD  Y  transmitter.  Solid  lines  are  actual  data,  circles  stand  for  NSMS  predictions . 70 

Figure  27:  Recovered  total  NSMS  from  calibration  BUD  measurements  for  M-75  (blue),  37-mm 

(green),  and  M-60  (red)  UXO . 70 

Figure  28.  Photo  and  schematic  diagram  of  the  MPV  sensor . 71 


Figure  29:  Multi-object  MPV  data  collection  setup  (right).  The  red  circle  corresponds  to  the  MPV 
head,  which  was  placed  stationary;  the  targets  were  moved  along  the  blue  line.  The  center  of 
the  first  target  (the  81 -mm)  was  placed  at  the  blue  points,  and  the  distance  between  the  first 
and  second  targets  was  kept  fixed . 72 

Figure  30:  Inverted  polarizability  principal  elements  for  two  targets  in  three  different  setups; 
results  for  the  81 -mm  projectile  at  left  and  for  the  40-mm  munition  at  right.  In  all  three  cases 
the  targets  were  horizontal,  and  the  vertical  distance  between  the  MPV  center  and  the  81 -mm 
was  40  cm.  The  center  to  the  center  coordinate  differences  between  the  8 1  -mm  and  40-mm 


projectiles  are  (-25,  0,  0)  cm,  (-40,  0,  0)  cm,  and  (-25,  0,  25)  cm . 73 

Figure  31:  Inverted  total  ONMS  for  81  mm  (left)  and  40  mm  (right)  projectiles  for  three  different 

cases . 73 

Figure  32:  Camp  Sibert  anomalies:  4.2  inch,  base  plates  and  partial  mortars . 75 

Figure  33:  Camp  Sibert  EM-63  near  field  distributions:  Left  and  middle  columns:  actual  and 

modeled  data  respectively.  Right  column:  misfits . 77 

Figure  34:  Inverted  total  NSMS  for  all  anomalies:  4.2"  mortars,  base  plates,  and  partial  mortars . 78 

Figure  35:  Left:  Classification  features.  Right:  ROC  curve  of  NSMS  performance . 78 

Figure  36:  a)  Unexploded  shell  from  Cell  No.  7  and  (b,c)  the  two  false  alarms  obtained  by  the 

SVM  classifier  using  A  and  Q(tl5)/  Q(t1 )  as  discriminators . 79 

Figure  37:  Result  of  the  SVM  classification  for  the  Camp  Sibert  anomalies  using  the  logarithms 


of  A  and  R  =  Q(tls)  /  Q(tt  ) .  The  SVM  has  been  trained  with  capacity  C  =  1 0  and  kernel  width  a 
=  1/200.  The  small  markers  denote  the  ground  truth  for  both  training  (hollow)  and  testing 


Sky  Research,  Inc. 


vi 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


(solid)  cells.  The  larger  markers  highlight  the  cases  where  there  is  disagreement  between  the 
ground  truth  and  the  SVM  prediction . 80 

Figure  38:  Result  of  the  SVM  classification  for  the  Camp  Sibert  anomalies  using  the  logarithms 
of  the  Pasion-Oldenburg  parameters  k  and  y.  The  SVM  here  has  a  capacity  C  =  9.  The  small 
markers  denote  the  ground  truth  for  both  training  (hollow)  and  testing  (solid)  cells.  The  larger 
markers  show  the  wrong  SVM  predictions . 82 

Figure  39:  Result  of  the  SVM  classification  for  the  Camp  Sibert  Anomalies  using  the  logarithms 
of  the  Pasion-  Oldenburg  parameters  J3  and  y.  The  SVM  capacity  C  =  105.  The  small  markers 
denote  the  ground  truth  for  both  training  (hollow)  and  testing  (solid)  cells.  The  larger  markers 
highlight  the  wrong  predictions  made  by  the  SVM . 82 

Figure  40:  SVM  classification  of  the  Camp  Sibert  Anomalies  using  the  logarithms  of  k,  b,  and  g. 

The  SVM  has  C  =  9.  The  small  markers  denote  the  ground  truth  for  both  training  and  testing 
cells.  The  larger  markers  highlight  the  cases  where  there  is  disagreement  between  the  ground 
truth  and  the  SVM  prediction . 83 

Figure  41:  Log-scale  plot  of  Q(tl5)  /  Qfy  vs.  k  for  Camp  Sibert  data  classification.  Left:  Ground 

truth.  Right:  X-means  clustering  result . 86 

Figure  42:  Classification  of  216  targets  into  five  classes  using  a  bivariate  normal  mixture.  Also 

shown  are  the  95%  confidence  ellipses . 86 

Figure  43:  Five  ROC  curves  that  indicate  the  performance  of  the  mixed  model  approach  to  Camp 

Sibert  data . 87 

Figure  44:  Found  Clutter  Items  on  SLO  UXO  live  sites . 88 

Figure  45:  Inverted  total  NSMS  time  decay  profiles  for  the  2.36"  partial  rocket.  The  green  lines 

depict  calibration  data  and  the  red  lines  correspond  to  blind  SLO  TEMTADS  data  sets . 90 

Figure  46:  Inverted  total  NSMS  time  decay  profiles  for  4.2"  mortars  (top  left),  81-mm  projectiles 
(top  right),  2.36"  rockets  (bottom  left),  and  60-mm  mortars  (bottom  right)  in  the  SLO 
TEMTADS  test.  The  green  lines  depict  calibration  data  and  the  red  lines  correspond  to  blind 
data  sets . 90 

Figure  47:  Result  of  the  supervised  clustering  classification  for  the  SLO-TEMTADS  anomalies 
using  the  logarithms  of  and  Maa{t&0) .  The  supervised  clustering  has  been 

trained  with  calibration  data.  The  red  markers  correspond  to  clutters  and  the  white  ones  to 

TOI . 91 

Figure  48:  ROC  curve  for  SLO  TEMTADS  test  data . 92 

Figure  49:  ROC  curve  for  SLO  MetalMapper  test  data . 93 

Figure  50:  ROC  for  SLO  TEMTADS  data  for  individual  TOI . 93 

Figure  51:  ROC  for  SLO  MetalMapper  data  sets:  individual  TOI . 94 

Figure  52:  60-mm  mortars  actually  found  in  calibration  cells  #410  and  #489 . 95 

Figure  53:  Left:  Principal  elements  of  the  polarizability  tensor  versus  time  for  a  60mm  mortar  in 
the  SLO  study.  Right:  Total  NSMS  time -decay  curves  for  the  same  cases.  The  red  curve 
corresponds  to  calibration  Cell  #489  and  the  blue  curve  to  calibration  Cell  #410 . 95 

Figure  54:  Comparison  between  library  and  inverted  blind  tests  for  the  dipole  model  (left)  and 

NSMS  model  (right) . 96 

vii 


Sky  Research,  Inc. 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Figure  55:  ROC  curves  for  SLO  TEMTADS  and  SLO  MetalMapper  discrimination  studies. 

Green  and  red  curves:  Sky/UBC  dipole  results;  blue  curve:  NSMS  results  obtained  by  our 
Dartmouth/Sky  group . 96 

Figure  56:  ROC  curves  for  SLO  BUD  discrimination  studies . 97 

Figure  57:  SLO  TEMTADS  test  Cell  #16.  Left:  All  25  eigenvalues  vs.  time.  Right:  Four  highest 
eigenvalues  vs.  time.  The  target  response  is  weak  and  mixed  with  the  sensor’s  electronic  and 
background  noise . 98 

Figure  58:  SLO  TEMTADS  test  Cell  #103.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above¬ 
threshold  eigenvalues  vs.  time.  Only  two  eigenvalues  are  above  the  threshold,  indicating  a 
low  signal-to-noise  ratio . 98 

Figure  59:  SLO  TEMTADS  test  Cell  #241.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above¬ 
threshold  eigenvalues  vs.  time.  There  more  than  three  eigenvalues  above  the  threshold,  which 
indicates  that  the  cell  contains  more  than  one  target.  The  curves  decay  fast,  illustrating  that 
the  targets  are  small . 99 

Figure  60:  SLO  TEMTADS  test  Cell  #441.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above¬ 
threshold  eigenvalues  vs.  time.  There  more  than  three  eigenvalues  above  the  threshold, 
indicating  that  the  cell  contained  more  than  one  target.  The  fast-decaying  curves  illustrate  that 
the  targets  have  thin  walls  or  are  small . 99 

Figure  61:  SLO  TEMTADS  test  Cell  #444.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above¬ 
threshold  eigenvalues  vs.  time.  There  more  than  three  eigenvalues  above  the  threshold, 
indicating  that  the  cell  contained  several  targets.  In  addition,  the  curves  decay  fast,  illustrating 
that  the  targets  are  small . 100 

Figure  62:  SLO  TEMTADS  test  Cell  #748.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above¬ 
threshold  eigenvalues  vs.  time.  More  than  three  fast-decaying  above-threshold  eigenvalues 
indicate  the  presence  of  several  small  targets . 100 

Figure  63:  SLO  TEMTADS  test  Cell  #1285.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above¬ 
threshold  eigenvalues  vs.  time.  Again,  the  eigenvalues  indicate  that  there  are  several  small 
targets  in  the  cell . 100 

Figure  64:  TEMTADS  multi-static  response  matrix  eigenvalues  versus  time  for  some  samples  of 

requested  anomalies . 102 

Figure  65:  TEMTADS  multi-static  response  matrix  eigenvalues  versus  time  for  a  105-mm  FIE 
projectile  and  a  105-mm  HEAT  round  (top  row),  an  M-48  Fuze  and  a  37-mm  munition 
(center  row),  and  two  clutter  scenarios,  one  with  two  items  (left)  and  another  with  several 
(right)  (third  row) . 103 

Figure  66:  Inverted  total  ONVMS  time -decay  profiles  for  four  Camp  Butner  targets:  (top  row) 

105-mm  HE  munition  and  105-mm  HEAT  round,  and  (bottom)  M-48  Fuze  and  37-mm 

projectile  with  copper  band . 105 

Figure  67:  Inverted  total  ONSMS  time  decay  profiles  for  a  37-mm  projectile  without  copper 

band . 105 

Figure  68:  ROC  curve  for  the  Camp  Butner  TEMTADS  test  data . 106 

Figure  69:  Left :  Scatter  plot  for  all  MM  anomalies  based  on  the  extracted  total  ONVMS.  Right : 

Probability  function  for  all  MM  anomalies . 107 

viii 


Sky  Research,  Inc. 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Figure  70:  Inverted  magnetic  dipole  polarizability  (left)  and  total  ONVMS  (right)  time-decay 
profiles  for  MM  anomaly  #2504.  The  thin  red  lines  show  a  library  sample,  while  the  thick 
blue  and  green  lines  show  the  inversion  results . 107 

Figure  7 1 :  Result  of  the  supervised  clustering  classification  for  the  Camp  Butner  MM  anomalies 
using  the  logarithms  of  and  M  (t ft  M  (t,0) .  The  supervised  clustering  was  trained 

with  calibration  data.  The  red  markers  correspond  to  clutter  and  the  green  ones  to  TOI . 108 

Figure  72:  ROC  curve  for  Camp  Butner  MetalMapper  test  data . 109 

Figure  73:  Left:  Total  ONVMS  time -decay  curves  for  a  105  mm  projectile  in  the  camp  Butner, 

NC  study.  Right:  Principal  elements  of  the  polarizability  tensor  versus  time  for  the  same 

case . 109 


Figure  74:  ROC  curve  for  Camp  Butner  live  site  classification.  100%  UXO  were  identified 
correctly,  with  only  295  false  positive  rate.  The  total  number  of  anomalies  is  2291.  The  blue 
dot  corresponds  to  a  threshold  in  the  dig  list,  when  the  boundary  between  UXO  and  clutter 
was  assumed  after  scoring.  The  cyan  dot  specifies  the  actual  position  of  this  boundary.  In 
ideal  circumstances  the  blue  and  cyan  points  will  coincide.  Performing  extra  digs,  however, 
helps  maintain  better  statistics  and  improve  the  results . 1 14 

Figure  75:  Camp  Butner  single -object  inverted  data  clustering.  Left:  Results  of  weighted-linkage 
clustering  using  Mahalanobis  distances  for  single-object  inverted  EMI  features.  Right:  All 
four  identified  UXO  (black)  after  a  second  clustering  within  a  smaller  domain 

(log#  e[2;8],  b  e[0.05;21,  g  e[0.05;2])  using  Ward  linkage  and  Euclidean  distances . 1 14 

Figure  76:  Camp  Butner  clusters  used  to  train  the  first  GM  classifier  and  its  results.  Left: 

Assumed  UXO  clusters  used  to  generate  the  3-component  GM  classifier.  Right:  Score 

histogram  showing  the  number  of  anomalies  scored  within  a  particular  range  of  the 

log(probability  density)  in  arbitrary  units.  The  ground  truth  was  requested  based  on 

thresholding  the  log(score)  at  the  externally  selected  value  of  ~0.5 . 1 16 

Figure  77:  Updated  GMM  classifier  after  confirming  118  UXO  in  the  Camp  Butner  data.  Left: 
GM-classifier  score  iso-surfaces  in  the  classification  case  based  on  all  currently  identified 
UXO  targets  (118  items),  in  the  feature  space  corresponding  to  3-object  EMI  inversion. 

Right:  Updated  score  histogram  showing  the  number  of  anomalies  scored  within  a  particular 
range  of  the  log(probability  density)  in  arbitrary  units.  An  additional  20  items  were  requested 
for  statistics  to  probe  the  region  corresponding  to  log(score)  within  [-6;  -5]  (this  region  was 
identified  with  external  input  from  an  expert  by  observing  the  corresponding  score  iso¬ 
surfaces  and  the  histogram  behavior) . 116 

Figure  78:  GMM  classifier  results  for  37-mm  targets.  Left:  1-component  GM-classifier  score  iso¬ 
surfaces  in  the  classification  case  based  solely  on  identified  37-mm  UXO  targets,  in  the 
feature  space  corresponding  to  3-object  EMI  inversion.  Right:  Score  histogram  showing  the 
number  of  anomalies  scored  within  a  particular  range  of  the  log(probability  density)  in 
arbitrary  units.  A  total  of  174  anomalies  were  requested  (with  118  being  already  known) 
based  on  the  log(score)  cut-off  value  of  about  -6  (specified  externally  by  an  expert  to  allow 
enough  statistics) . 116 

Figure  79:  GMM  classifier  for  48-mm  targets.  Left:  1-component  GM-classifier  score  iso-surfaces 
in  the  classification  case  based  solely  on  identified  48-mm  UXO  targets,  in  the  feature  space 
corresponding  to  3 -object  EMI  inversion.  Right:  Score  histogram  showing  the  number  of 
anomalies  scored  within  a  particular  range  of  the  log(probability  density)  in  arbitrary  units.  A 
total  of  53  anomalies  were  requested  (with  27  being  already  known)  based  on  the  log(score) 
cut-off  value  of  about  -5  (specified  externally  by  an  expert  to  allow  enough  statistics) . 117 


Sky  Research,  Inc. 


IX 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Figure  80:  GMM  classifier  for  105-mm  targets.  Left.  1-component  GM-classifier  score  iso¬ 
surfaces  in  the  classification  case  based  solely  on  identified  105-mm  UXO  targets,  in  the 
feature  space  corresponding  to  3-object  EMI  inversion.  Right  Score  histogram  showing  the 
number  of  anomalies  scored  within  a  particular  range  of  the  log(probability  density)  in 
arbitrary  units.  A  total  of  36  anomalies  were  requested  (with  18  being  already  known)  based 
on  the  log( score)  cut-off  value  was  about  -20  (specified  externally  by  an  expert  to  allow 
enough  statistics) . 1 17 

Figure  81:  Final  GMM  classifier  for  Camp  Butner.  Left:  3-component  GM-classifier  score  iso¬ 
surfaces  in  the  classification  case  based  on  all  identified  UXO  targets,  in  the  feature  space 
corresponding  to  3 -object  EMI  inversion.  Right:  Score  histogram  showing  the  number  of 
anomalies  scored  within  a  particular  range  of  the  log(probability  density)  in  arbitrary  units.  A 
total  of  377  anomalies  were  scored  as  UXO  based  on  the  log(score)  cut-off  value  of  about  - 
10  (specified  externally  by  an  expert  to  allow  enough  statistics  and  iso-surface  separation 
from  identified  UXO  clusters) . 117 


Sky  Research,  Inc. 


x 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


List  of  Tables 

Table  1 :  MetalMapper  receiver  locations  with  respect  to  the  center  of  the  Z  transmitter  loop . 55 

Table  2:  Inverted  location  and  orientation  for  TEMTADS  data . 61 

Table  3:  BUD  receiver  locations  with  respect  to  the  origin . 68 

Table  4:  SVM  classification  of  Camp  Sibert  anomalies  using  k  and  R  with  C  =  10 . 80 

Table  5:  SVM  classification  of  Camp  Sibert  anomalies  using  y  and  k  with  C  =  9 . 81 

Table  6:  SVM  classification  of  Camp  Sibert  anomalies  using  p  and  y  with  C  =  105 . 83 

Table  7:  SVM  classification  of  Camp  Sibert  anomalies  using  the  complete  NSMS  time  decay . 83 


Sky  Research,  Inc.  xi  January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


List  of  Acronyms 

3D  Three-Dimensional 

AIC  Akaike  Information  Criterion 

APG  Aberdeen  Proving  Ground 

BIC  Bayesian  Information  Criterion 

BOR  Body  of  Revolution 

BUD  Berkeley  UXO  Discriminator 

cm  Centimeter 

CRREL  Cold  Regions  Research  and  Engineering  Laboratory 

DE  Differential  Evolution 

DGPS  Differential  Global  Positioning  System 

DoD  Department  of  Defense 

DOE  Department  of  Energy 

EM  Expectation  Maximization 

EMI  Electromagnetic  Induction 

ESTCP  Environmental  Security  Technology  Certification  Program 

FD  Frequency  Domain 

GSEA  Generalized  Standardized  Excitation  Approach 

HAP  (1)  magnetic  field  vector  H,  (2)  vector  potential  A,  (3)  scalar  magnetic  Potential  y/ 

Hz  Hertz 

ID  Identification 

LS  Least  Squares 

m  Meter 

MAS  Method  of  Auxiliary  Sources 

MD  Multi-Dipole 

xii 


Sky  Research,  Inc. 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


MDL 

MEG 

ML 

mm 

ms 

Ms 

MM 

MNM 

MPV 

MT 

MUSIC 

NMMS 

NRL 

NSMC 

NSMS 

ONVMS 

PCA 

PNN 

RBF 

ROC 

RSS 

RTS 

SEA 

SFS 

SLO 

SNR 


Minimum  Description  Length 

Magnetoencephalographic 

Maximum  Likelihood 

Millimeter 

Millisecond 

Microsecond 

MetalMapper 

Multivariate  Normal  Mixture 
Man-Portable  Vector 
Mixed  Theory 

Multiple  Signal  Classification 

Non-metric  Multidimensional  Scaling 

Naval  Research  Laboratory 

Normalized  Surface  Magnetic  Charge 

Normalized  Surface  Magnetic  Source 

Ortho-Normalized  Volume  Magnetic  Source  Model 

Principal  Component  Analysis 

Probabilistic  Neural  Network 

Radial  Basis  Function 

Receiver  Operating  Characteristic 

Residual  Sum  of  Squares 

Robotic  Total  Station 

Standard  Excitation  Approach 

Scattered  Field  Singularities 

San  Luis  Obispo 

Signal-to-Noise  Ratio 


xiii 


Sky  Research,  Inc. 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


ss 

Sum  of  Squares 

SVM 

Support  Vector  Machine 

TD 

Time  Domain 

TEM 

Time  Domain  Electromagnetic 

TEMTADS 

Time  Domain  Electromagnetic  Towed  Array  Detection  System 

TNLL 

Twice-Negative  Log-Likelihood  Function 

TOI 

Target  Of  Interest 

TSA 

Thin  Skin  Approximation 

UXO 

Unexploded  Ordnance 

VRM 

Viscous  Magnetic  Remanence 

Sky  Research,  Inc, 


xiv 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Abstract 

The  research  described  in  this  report  was  conducted  in  fulfillment  of  Project  MM- 1572,  “A 
Complex  Approach  to  UXO  Discrimination:  Combining  Advanced  EMI  Forward  and  Statistical  Signal 
Processing,”  submitted  to  the  Strategic  Environmental  Research  and  Development  Program  (SERDP)  in 
response  to  the  Munitions  Management  Statement  of  Need  MMSON-07-04  “Advanced  Technologies  for 
Detection,  Discrimination,  and  Remediation  of  Munitions  and  Explosives  of  Concern  (MEC):  UXO 
Technology”. 

The  primary  objective  of  the  project  MM- 1572  was  to  develop  and  validate  innovative,  robust, 
and  practical  approaches  for  UXO  localization  and  classification  under  realistic  (noisy,  cluttered 
background)  field  conditions  by  combining  advanced  electromagnetic  induction  (EMI)  forward  and 
statistical  signal  processing  methodologies.  In  a  real  field  the  electromagnetic  signals  become  convoluted 
with  noise  due  to  the  instrument,  magnetic  soil  and  widespread  background  clutter.  Therefore,  the 
rationale  for  a  statistical  approach  is  to  use  an  advanced  statistical  approach  to  reduce  the  impact  of  these 
noises  to  minimum.  This  project  provides  mathematical  fundamentals,  physical  meanings  and  practical 
realizations  of  forward,  inverse  and  statistical  signal  processing  approaches  for  unexploded  ordnance 
(UXO)  detection  and  discrimination  at  live-UXO  sites.  Namely,  under  this  project  first  we  developed  and 
implemented  advanced,  physically  complete  forward  EMI  models  such  as,  the  normalized  surface 
magnetic  source  (charge/dipole)  model  (NSMS),  and  ortho-normalized  volume  magnetic  source 
(ONVMS)  technique  for  accurately  representing  the  EMI  responses  of  subsurface  metallic  targets,  then 
we  combined  these  advanced  models  with  EMI  data  inversion  approaches,  such  as  the  gradient  search, 
direct  search-differential  evolution  and  etc,  for  advanced  EMI  sensor  data  inversion;  third  we  extended 
the  advanced  statistical  signal  processing  approaches,  i.e.  support  vector  machines,  Gaussian  mixture 
models,  for  discriminating  UXO  targets  from  non-hazardous  anomalies.  Finally,  the  combined  advanced 
EMI  forward  and  statistical  models  were  applied  to  ESTCP  live  site  UXO  data  sets.  Live  site 
discrimination  studies  showed  the  excellent  discrimination  performance  of  the  advanced  models  when 
applied  to  next-generation-sensor  data  collected  at  various  live  sites,  such  as  Camp  Sibert,  AL,  San  Luis 
Obispo  (SLO),  CA,  and  Camp  Butner,  NC  as  well  as  APG  test  sites.  The  technology  was  able  to  single 
out  UXO  ranging  in  caliber  from  25  mm  up  to  155  mm.  In  addition,  the  ONVMS  technique  was  seen  to 
provide  excellent  classification  in  both  single-  and  multiple-target  scenarios  when  combined  with 
advanced  multi-axis/transmitter/receiver  sensors  data. 
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1  Introduction 

1.1  Background  and  objectives 

The  research  described  in  this  report  was  conducted  in  fulfillment  of  Project  MM- 1572,  “A 
Complex  Approach  to  UXO  Discrimination:  Combining  Advanced  EMI  Forward  and  Statistical  Signal 
Processing,”  submitted  to  the  Strategic  Environmental  Research  and  Development  Program  (SERDP)  in 
response  to  the  Munitions  Management  Statement  of  Need  MMSON-07-04,  “Advanced  Technologies  for 
Detection,  Discrimination,  and  Remediation  of  Munitions  and  Explosives  of  Concern  (MEC):  UXO 
Technology.” 

The  well-known  and  prohibitive  cost  of  carefully  excavating  all  geophysical  anomalies  detected 
at  lands  contaminated  with  unexploded  ordnance  (UXO)  is  one  of  the  greatest  impediments  to  performing 
an  efficient  and  thorough  cleanup  of  former  battlefields  and  of  Department  of  Defense  (DoD)  and 
Department  of  Energy  (DOE)  sites.  Innovative  discrimination  techniques  are  required  that  can  quickly 
and  reliably  distinguish  between  hazardous  UXO  and  non-hazardous  metallic  items.  The  key  to  success 
lies  in  the  development  of  advanced  processing  techniques  that  can  analyze  and  process  sophisticated 
magnetic  or  electromagnetic  induction  data,  with  its  novel  waveforms,  ever  improving  quality,  and  vector 
or  tensor  character,  so  as  to  maximize  the  probability  of  correct  classification  and  minimize  the  false- 
alarm  rate. 

The  main  objective  of  the  proposed  work  was  to  combine  physically  complete  forward  models 
and  state-of-the-art  statistical  signal  processing  methodologies  to  carry  out  dependable  and  robust  UXO 
discrimination  in  difficult  and  noisy  sites  starting  from  data  provided  by  current  electromagnetic 
induction  (EMI)  sensors.  Specifically,  our  objectives  were  to 

•  Extend  physically  complete  forward  approaches  like  the  normalized  surface  magnetic  source 
model  and  the  orthonormalized  volume  magnetic  source  model  to  treat  realistic  data  sets. 

•  Formulate  and  develop  an  inversion  framework  featuring  robust  regularization  and  parameter- 
determination  methodologies  (for  both  linear  intrinsic  signatures  and  non-linear  extrinsic 
particulars)  based  on  advanced  signal  processing  algorithms. 

•  Combine  EMI  models  and  statistical  methodologies  to  process  complex,  heterogeneous 
geophysical  data. 

•  Demonstrate  the  discrimination  capability  of  the  combined  approach  by  applying  it  to  blind  live- 
site  UXO  discrimination  studies. 
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1.2  Report  structure 

Chapter  2  outlines  the  theoretical  basis  of  the  advanced  EMI  forward  models  that  we  use  to 
represent  the  EMI  response  of  obscured  targets,  the  workhorses  at  the  core  of  our  whole  procedure.  We 
first  present  the  single-dipole  model,  which  is  usually  insufficient  in  itself,  but  undergirds  the  others.  We 
then  introduce  and  study  in  detail  the  NSMS  model,  which  distributes  dipoles  on  a  closed  surface 
surrounding  a  target  of  interest.  After  that  we  derive  and  describe  the  ONVMS  technique,  which  infuses 
dipoles  throughout  the  subsurface  volume  illuminated  by  a  sensor.  We  end  by  describing  a  data- 
preprocessing  technique  based  on  joint  diagonalization  that  estimates  the  number  of  targets  in  a 
measurement  with  no  need  for  data  inversion;  the  method,  moreover,  can  provide  initial  estimates  of 
tai'gct  locations  and  perform  rudimentary  discrimination. 

Chapter  3  discusses  inverse  models:  the  methods  used  to  harness  the  forward  models  so  they 
provide  relevant  intrinsic  and  extrinsic  information  starting  from  measured  data.  After  presenting  some 
traditional  gradient-search  based  methods  and  pointing  out  some  of  their  limitations  we  describe 
differential  evolution,  a  state-of-the-art  global-search  method,  similar  in  character  to  genetic  algorithms, 
that  has  shown  remarkable  flexibility  and  usefulness.  We  end  by  describing  the  HAP  method,  a  semi- 
analytic  non-iterative  procedure  to  locate  buried  targets. 

Chapter  4  introduces  the  theory  and  application  of  various  statistical  techniques  for  UXO 
classification.  First,  the  Mixed  Models  for  UXO  discrimination  is  presented  and  the  discrimination  based 
on  NSMS  data  in  discussed.  Finally,  some  other  supervised  (SVM)  and  unsupervised  classification 
techniques  are  discussed  in  application  for  UXO  discrimination. 

Chapter  5  presents  and  describes  the  next-generation  EMI  sensors — the  MetalMapper,  the 
TEMTADS  array,  the  MPV  portable  instrument,  and  the  BUD  system — that  took  all  the  data  we  use  and 
that  represent  the  state  of  the  art  in  UXO  remediation  hardware.  We  present  the  results  of  several  testing 
and  validation  studies  carried  out  on  laboratory,  test  stand  and  US  army  standardized  Aberdeen  Proving 
Ground  in  Maryland  test-site  data  from  these  devices. 

Finally,  Chapter  6  provides  a  detailed  account  of  the  discrimination  and  classification  studies 
performed  on  data  from  actual  UXO  sites — the  Camp  Sibert  in  Alabama,  Camp  San  Luis  Obispo  in 
California,  and  Camp  Butner  in  North  Carolina — in  which  several  combinations  of  the  techniques 
presented  in  the  previous  chapters  were  used.  We  describe  our  solution  strategies  and  the  results  we 
obtained. 
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2  Forward  models 

2.1  Introduction 

UXO  discrimination  is  an  inverse  problem  that  demands  a  fast  and  accurate  representation  of  a 
target’s  EMI  response.  Much  of  our  research  in  this  project  has  had  to  do  with  the  development, 
implementation,  and  honing  of  models  that  provide  such  representations  in  a  physically  complete,  noise- 
tolerant  way  that  allows  them  to  perform  adequately  in  realistic  settings  and  to  set  the  stage  for 
dependable  live-site  UXO  discrimination.  In  this  section  we  provide  an  overview  of  the  theoretical  basis 
and  implementation  details  of  several  models  of  increasing  usefulness  and  sophistication,  starting  with  the 
standard  point-dipole  approach,  continuing  with  the  generalized  standardized  excitation  approach  (GSEA) 
and  the  normalized  surface  magnetic  source  (NSMS)  model,  and  culminating  with  the  orthonormalized 
volume  magnetic  source  (ONVMS)  model.  We  finish  the  section  by  describing  a  fast,  inversion-free 
method  to  estimate  the  number  of  buried  targets. 

The  most  frequently  used  method  for  representing  the  EMI  response  of  a  metallic  target  in  both 
frequency  and  time  domains  approximates  the  whole  object  with  a  set  of  orthogonal  co-located  point 
dipoles  that  fire  up  in  response  to  the  primary  field;  the  induced  dipole  moment  is  related  to  the  primary 
field  through  a  symmetric  polarizability  tensor.  The  use  of  this  dipole  approximation  is  motivated  by  its 
speed  and  simplicity;  this  simplicity,  however,  rests  on  assumptions  that  often  become  problematic  and 
limit  the  model’s  usefulness.  One  such  assumption  is  that  the  buried  target  of  interest  is  either  far  enough 
from  the  transmitter  loop,  or  small  enough,  that  the  primary  field  is  essentially  uniform  throughout  its 
extent.  Usually,  complex  targets  composed  of  different  materials  and  different  sections  that  contribute 
appreciably  to  the  response — and,  in  the  case  of  UXO,  containing  such  complicating  features  as  fins  and 
rings — simply  cannot  be  modeled  accurately  with  a  single  point  dipole.  Such  cases  require  more 
advanced  methods  that  will  capture  the  underlying  physics  correctly.  One  such  technique  is  the  NSMS 
model. 

The  NSMS  method  [1-4]  can  be  considered  as  a  generalized  surface  dipole  model,  and  indeed 
reduces  to  the  point  dipole  model  in  a  special  limiting  case.  The  NSMS  approach  models  an  object’s 
response  to  the  primary  field  of  a  sensor  by  distributing  a  set  of  equivalent  elementary  magnetic 
sources — normally  oriented  dipoles  in  this  case — over  an  auxiliary  surface  that  surrounds  it.  Such  a 
surface  distribution  can  be  hypothetically  generated  by  spreading  positive  magnetic  charge  over  the  outer 
side  of  the  equivalent  surface  (usually  a  prolate  spheroid)  and  an  identical  distribution  of  opposite  sign  on 
its  inner  side  [5],  resulting  in  a  double  layer  of  magnetic  charge  separated  by  an  infinitesimal  distance. 
This  double  layer  introduces  the  proper  discontinuities  in  the  tangential  components  of  the  magnetic  flux 
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density  vector  B  but  does  not  affect  the  transition  of  its  normal  component,  which  must  always  be 
continuous  given  the  lack  of  free  magnetic  charges  in  nature.  The  resulting  magnetic -moment  distribution 
radiates  a  field  that  by  construction  satisfies  the  governing  EMI  equations  and  can  thus  account  for  the 
secondary  field  outside  the  object.  The  particulars  of  location  and  orientation  are  divided  out  by 
normalizing  the  dipole  density  at  every  point  with  the  component  of  the  primary  magnetic  field  normal  to 
the  surface.  The  resulting  surface  amplitude  Q  of  the  NSMS  distribution  is  a  property  of  the  object,  and 
its  integral  Q  over  the  surface  constitutes  a  sort  of  global  magnetic  polarizability  that  is  independent  of 
the  computational  constructs — primary  field,  surrounding  surface,  object  location  and  orientation,  etc. — 
introduced  for  its  determination.  The  surface  amplitude  can  be  determined  directly  for  library-matching 
puiposes  by  minimizing  the  difference  between  measured  and  modeled  data  for  a  known  combination  of 
object  and  sensor  at  a  given  relative  location  and  orientation. 

The  NSMS  technique  has  demonstrated  good  computational  speed  and  superior  classification 
performance  when  applied  to  EMI  datasets  consisting  of  well-isolated  single  targets,  but  is  found  to 
degrade  quickly  on  both  counts  when  confronted  with  multi-target  cases.  This  has  forced  us  to  generalize 
the  model  further  and  develop  the  ONVMS  procedure. 

The  ONVMS  model,  a  further  extension  of  NSMS,  is  based  on  the  assumption  that  a  collection  of 
subsurface  objects  can  be  replaced  with  a  set  of  magnetic  dipole  sources,  distributed  over  a  volume. 
Since  all  actual  radiating  sources  are  located  within  the  scatterers — rather  than  in  the  soil  or  air — the 
spatial  distribution  of  these  fictitious  dipoles  (their  amplitudes  scaled  by  the  primary  field)  indicates  the 
locations  and  orientations  of  any  targets  present  inside  the  computational  volume.  The  great  advantage  of 
the  ONVMS  technique  over  the  other  models  discussed  above  is  that  it  takes  into  account  mutual 
couplings  between  different  sections  of  the  different  targets  while  simultaneously  avoiding  the  appearance 
of  singular  matrices  in  multi-target  situations.  It  is  thus  gracefully  indifferent  to  the  number  of  targets: 
Once  the  amplitudes  and  the  locations  of  the  corresponding  dipoles  are  determined,  one  need  only  look  at 
their  clustering  patterns,  compute  the  time -dependent  total  polarizability  tensor  for  each  group,  and 
subsequently  diagonalize  each  such  tensor  using  joint  diagonalization  (see  immediately  below  for  another 
application).  The  resulting  diagonal  elements  have  been  found  to  be  intrinsic  to  the  objects  they  represent, 
and  can  be  used,  on  their  own  or  combined  with  other  quantities,  in  discrimination  analysis.  Recent 
ESTCP  live-site  discrimination  studies  have  clearly  indicated  the  superior  discrimination  performance 
(illustrated  in  chapter  6)  of  the  ONVMS  method  in  combination  with  the  statistical  processing  approaches 
described  in  Chapter  4. 

One  of  the  main  challenges  one  faces  when  attempting  multi-target  inversion  and  classification  is 
the  inability  to  estimate  the  number  of  targets.  Under  this  project  we  implemented  a  technique  based  on 
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joint  diagonalization  that  estimates  the  number  of  targets  present  in  the  field  of  view  of  the  sensor  as  it 
takes  a  data  shot,  in  real  time  and  without  requiring  a  forward  model,  and,  in  a  good  number  of  cases, 
even  provides  the  capability  to  perform  a  quick  inversion-free  characterization  and  classification  of  these 
targets.  JD  determines  the  eigenvalues  and  eigenvectors  of  a  square  time-  or  frequency-dependent  multi¬ 
static  response  (MSR)  matrix  synthesized  directly  from  measured  data.  The  number  of  nonzero 
eigenvalues  of  the  matrix  (i.e.,  those  above  a  noise  threshold)  is  related  to  the  number  of  elementary 
sources  in  the  illuminated  cell;  moreover,  the  time-decay  patterns  of  these  non-vanishing  eigenvalues  are 
intrinsic  properties  of  the  targets  to  which  the  sources  correspond  and  can  ultimately  provide  dependable 
classification  features. 

2.2  The  single -dipole  approximation 

According  to  the  Huygens  Equivalence  Principle,  an  object’s  entire  response  to  a  given  excitation 
can  be  approximated  as  the  summation  of  magnetic  fields  produced  by  elementary  magnetic 
dipoles/charges  placed  on  a  closed  surface  surrounding  the  target.  Using  the  superposition  principle,  this 
set  of  dipoles  can  be  approximated  as  one  independent  aggregate  dipole.  In  the  simple  dipole  model,  the 
secondary  magnetic  field  at  r  due  to  a  dipole  of  moment  m  is: 

H  =  — !—  (3RR-T)  m  =  5  m  (1) 

4  ttR3 

where  R  is  the  unit  vector  along  R  =  r'-rrf,  r  is  the  dipole’s  position,  and  I  is  the  identity  dyad  (see 
Figure  1).  The  dipole  moment  m  induced  by  the  primary  magnetic  field  Hpr  is  given  by 

m  =  I^[  Hpr(r',rd),  (2) 

where  M,  the  target's  magnetic  polarizability  tensor,  is  a  symmetric  matrix:  M  =M  ,  a,j3  =  x,y,z  . 
This  tensor  depends  on  the  scatterer’s  shape,  size,  and  material  properties.  In  a  coordinate  system  aligned 
with  the  scatterer’s  principal  axes  for  different  primary  magnetic  fields  Hpr(r,r ),  (2)  can  be  written  in 
matrix  form  as 

[m]  =  I^{Hpr].  (3) 
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Figure  1:  A  dipole’s  location  in  a  global  coordinate  system. 


Thus  the  secondary  magnetic  field  is 


H  =  G  - M-[hp‘  ]  =  [T] -[m]7"  , 


(4) 


where  [M]  is  a  1  x  6  dimensional  vector  whose  components  (Mxx,  Mxy,  Mxz,  Myy,  Myz,  MZ7)  correspond  to 
the  elements  of  the  target’s  magnetic  polarizability  tensor  M  and  [T]  is  a  3  x  6  matrix, 
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whose  elements  are  as  follows: 
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Once  the  vector  M  is  determined  the  magnetic  polarizability  tensor  M  is  constructed  as 
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and  finally  the  M  tensor’s  principal  polarizability  elements  are  determined  in  the  target  frame  coordinate 
system,  which  is  related  to  the  global  coordinate  system  via  the  Euler  rotation  tensor  A( i//,  0, (j>) ,  as 


'  Pr.  0 
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Body-of-revolution  (BOR)  symmetry  (which  most  UXO  possess)  dictates  that  Px  =  p  and  that  the  third 
Euler  angle  i//  is  zero.  We  thus  obtain 


A  = 


cos  0  cos  (j)  cost?  sin  ^  -sint? 

-sin^  cos  (f>  0 

sin^cos^  sin(?sin^  cos^ 


(8) 


where  6  and  </>  are  the  angles  between  the  local  and  global  axes.  Note  that  the  tensor  M  depends  on 
time  or  frequency  while  the  Euler  tensor  does  not.  This  suggests  that  one  could  apply  joint 
diagonalization  (along  the  lines  of  the  procedure  introduced  in  Section  2.5)  to  separate  the  polarizability 
eigenvalues  from  the  rotational  eigenvectors;  the  attitude  angles  can  in  turn  be  extracted  from  the  latter. 


2.3  Normalized  surface  magnetic  source  model 

2.3.1  Governing  equations 

Outside  the  objects  of  interest,  both  the  primary  and  the  secondary  magnetic  fields  in  air  and 
ground  (Hpr  and  Hsc  respectively,  both  measured  in  A/m)  are  irrotational  and  can  be  represented  by  scalar 
magnetic  potentials  (with  dimensions  of  current)  through 
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Hout  =  Hpr  +  Hsc  -  -V y/oui  =  -V^pr  -  V^sc .  (9) 

Inside  a  metallic  target  the  EMI  field  is  governed  by  a  diffusion  equation  and  can  be  expressed  as 


Hm  =  — VxAm, 

where  the  vector  potential  A  (also  with  dimensions  of  current)  obeys 


(10) 
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and  /j"1  —  ju'" ju0 ;  //„  —  An  -\()  [H/mJ  and  d"  are  respectively  the  magnetic  permeability  and  the  electric 
conductivity  of  the  object. 

The  total  interior  and  exterior  magnetic  fields  are  connected  at  the  surface  of  an  object  by  the 
customary  boundary  conditions 


hx(Hpr  +  Hsc)  =  rixHin, 
n  •  jU°ut  (Hpr  +  Hsc )  =  n  •  /t‘nHm , 


(12) 


where  h  is  a  unit  vector  normal  to  the  surface  and  //"  and  //"L"  are  the  relative  magnetic  permeabilities 

of  the  target  and  the  surrounding  medium  respectively.  Dividing  the  scatterer’s  surface  into  subsurfaces 
(patches)  allows  us  to  rewrite  the  boundary  conditions  in  the  following  convenient  form: 

[G][  P]  =  [V],  (13) 


where 


GP"  n  GP“  n  GP' 

n  r  r  n  r  r  n 

p 

n 

~Hpr' 

n 

GP  GP  GP 

u  u  u 

,  [P]= 

P 

u 

,  and  [V]  =  - 

Hpr 

u 

1 

1 _ 

P 

V 

Hw 

V 

(14) 


The  vector  [P]  contains  the  amplitudes  of  induced  magnetic  dipoles  oriented  along  the  n,  u  and  v 
directions,  where  n  is  normal  to  the  surface  and  v  =  n  x  u  and  u  =  v  x  n  are  tangential  to  it.  The 
impedance  matrix  [G]  depends  only  on  the  target’s  geometry  and  electromagnetic  parameters  and  not  on 
its  location  and  orientation;  its  elements  represent  the  exterior  and  interior  solutions,  expressed  in  terms  of 
dipole  sources  distributed  over  the  surface  using  a  Green  function  of  the  form  ejkR / AnR .  More  explicit 
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P  S'  /V  /V  /V 

forms  of  the  matrices,  where  2j  =  n,  u,  v,  are  presented  in  [4].  Finally,  [V]  is  the  primary  magnetic 

field  that  can  be  considered  as  a  driving/forcing  vector  for  the  [P] .  Since  most  actual  UXO  are 
heterogeneous  objects,  numerical  methods  may  require  tens  or  even  hundreds  of  thousands  of  unknowns 
for  magneto  quasi  static  boundary-value  problems.  This  makes  rigorous  3D-EM1  numerical  models 
impractical  for  UXO  discrimination.  To  achieve  the  goal  of  reproducing  EMI  responses  of  realistic 
objects  quickly,  accurately,  and  with  only  a  few  model  parameters  we  developed  the  normalized  surface 
magnetic  source  model. 


2.3.2  Theoretical  basis  of  NSMS 

The  NSMS  model  is  based  upon  the  assumption  that  the  entire  scatterer  can  be  replaced  with  an 
auxiliary  very  thin  surface  shell.  The  primary  magnetic  field  strikes  the  shell  and  induces  on  it  a  surface 
magnetization,  in  terms  of  which  the  secondary  scalar  potential  can  be  written  as  [5] 


^sc(r)  =  —  I  M(r')-V’-ds'  (15) 

4 n®s  R 

Here  R  =  r  -  r' ,  where  r  is  the  observation  point  and  r'  is  on  the  surface  S,  and  M(r')  is  a  surface 
density  of  magnetization,  which  can  be  defined  as  the  induced  magnetic  moment  per  unit  surface: 

m  =  .  The  surface  density  M  of  magnetic  polarization  may  be  resolved  at  every  point  on  S 

into  normal  and  tangential  components  by  means  of  the  identity 

M  =  (n  •  M)n  +  (n  x  M)  x  n ,  (16) 


and  combining  (16),  (15)  and  (9)  we  get  for  the  total  scattered  magnetic  field 

Hsc(r)  =  — —  V(6  (h'  ■  M(r ') )n'  -V'  —  ds' 

An  J  s  R 

— —  Vc6  (n'  x  M(r'))  xn'-V'-  ds'. 
An  rs  R 


(17) 


The  first  integral  in  (17)  may  be  interpreted  as  a  scalar  potential  due  to  a  double  layer  of  moment 

x(r')  =  (n'  •  M(r'))n'  =  crHI  (r')n' ,  (18) 

and  the  second  may  be  interpreted  as  a  scalar  potential  due  to  a  “free”  magnetic  charge  distribution 
proportional  to  a  discontinuity  in  the  normal  components  of  magnetic  flux.  Since  the  normal  component 
of  the  magnetic  field  is  always  continuous  across  a  boundary  between  two  media,  the  total  scattered 
magnetic  field  can  thus  be  written  as 
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Hsc(r)  =  -V4^,(r')g(r,r')^'  (19) 

where 

g(r,r')  =  -^-n'-V'  — l—  (20) 

An  Ir-r  I 


Thus  the  EMI  response  of  a  permeable  and  conducting  metallic  object  can  be  represented  using  a 
surface  density  crm(s').  At  every  point,  the  magnetic  flux  density  B  is 


B  =  /yH  +  M). 


(21) 


Using  Gauss’s  law  for  the  magnetic  flux  density  in  the  volume  enclosed  by  S  and  using  the  divergence 
theorem  we  obtain 


J  ' V  •  Bdv  =  ju0  (H  •  n'  +  M  •  n')ds' 

=  =  0, 


(22) 


and  it  follows  that  the  magnetization  density  at  a  given  point  on  the  surface  equals 


=  -//pV)(l+/V)), 


(23) 


where  Pis ')  is  in  general  position-dependent  on  S  surface.  In  other  words,  the  surface  magnetic  charge  is 
proportional  to  the  normal  component  H^(s’)  of  the  primary  magnetic  field.  This  motivates  us  to 
introduce  a  normalized  surface  distribution  Q(s ')  through 

a;>>-QG')[HpV)-n'],  (24) 


which  would  result  from  exciting  each  patch  of  the  surface  S  with  a  nonphysical  unit  primary  magnetic 
field  in  the  normal  direction.  After  combining  (19)  and  (24),  the  total  scattered  magnetic  field  can  be 
expressed  as 


Hsc  (r)  =  (j)s  Q(r')  [Hpr(r')  •  n']  Vg(r,r')ds' 


=  cj^Q(r')  [Hpr(r')-h'] 


3R(R  n'l-fl  n' 
4nR5 


ds'. 


(25) 


In  the  following  we  will  argue  that  Q.  and  in  particular  its  integral  over  the  surface, 

Q  =  ^nds' , 


(26) 
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contains  all  the  information  about  an  object  that  could  be  of  need  in  the  UXO  discrimination  problem, 
incorporating  the  effects  of  heterogeneity,  interaction  with  other  objects,  and  near-  and  far-field  effects. 
We  note  that  Q  has  dimensions  of  volume,  which  makes  it  comparable  to  the  polarizability  tensor 
elements  of  the  point  dipole  model  [6-13], 

2.3.3  Formulation  for  bodies  of  revolution;  determining  NSMS  amplitudes  from  data 

Most  UXO  are  bodies  of  revolution  (BOR),  and  the  simplicity  and  efficiency  afforded  by  this 
simplification  motivates  specializing  the  above  analysis  to  scatterers  with  BOR  symmetry.  The  best 
choice  for  auxiliary  surface  is  a  prolate  spheroid,  since  it  has  BOR  symmetry  but  at  the  same  time  has  the 
elongated  shape  of  UXO  and  can  be  made  to  have  a  definite  orientation.  We  take  a  spheroid  of  semiminor 
and  semimajor  axes  a  and  b  =  ea  with  e  >  1.  In  the  prolate  spheroidal  coordinate  system  (2,  //,  (p)  we  can 
write  (25)  in  the  form  (see  Figure  2) 

Hsc (r)  =  \  t hjdrj'l  ht,d<p' 0(77', 4 , <p') [H|ir (r') •  £] Vg(r, r ') ,  (27) 

where  the  prolate  spheroidal  coordinates  obey  -1  <  7  <  1 ,  0  <  2  < 00  >  0  <  <  2tt  ,  r  is  the  observation 

point,  h  fj  and  h,„  are  the  metric  coefficients 

\  =  f  *»  =  f  VO-^XC-D,  (28) 

the  spheroid  is  characterized  by  4  =  e  /  V e1  - 1  ,  and  d  —  2\jb2  -  a2  is  the  focal  distance. 

For  a  body  with  BOR  symmetry  the  NSMS  amplitude  is  azimuthally  constant,  and  moreover  the 
variation  of  the  induced  magnetic  charge  density  a  is  accounted  for  by  the  normal  component  of  the 
primary  magnetic  field.  This  implies  that  Q.(r/' ,%Q,cp')  =  Cl(r/') .  For  convenience  we  define 

Hsc(r)  =  D.(rj'yK(r/' ,r)drj' ,  (29) 

where 

K(/7',r)  =  jn2T[Hpr(r')  •  |']  Vgir^Ti'^h^dcp' ,  (30) 

and  assume  that  the  NSMS  can  be  approximated  by  a  series  of  expansion  functions  FJjf )  such  that 

M 

a(,')  =  Z£2.F(,')  (3D 

m= 1 
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Figure  2:  The  NSMC  that  are  distributed  on  a  prolate  spheroidal  surface  is  implemented  for  a  body  of  revolution. 
The  prolate  spheroidal  coordinate  system  is  specified  by  (£,  77,  <p). 


For  computational  simplicity,  in  the  subsequent  analysis  we  assume  the  expansion  functions  Fm(  77)  are  a 
set  of  orthogonal  pulse  functions  given  by 

1,  77'  GA77 

Fm{? 7')=  m  (32) 

0,  otherwise. 


The  expansion  in  terms  of  pulse  functions  is  a  “stairstep”  approximation  to  the  NSMS  distribution  on  the 
spheroid  along  77',  where  the  spheroidal  surface  is  divided  into  M  belts.  The  expansion  coefficient  Q„,  thus 
corresponds  to  the  NSMS  amplitude  at  the  777-th  belt.  Substituting  into  (29)  we  obtain 

1  M 

Hsc(r)  =  f  _Xn«Fm(ri')K(rf,T)drf ,  (33) 

m=  1 

and  the  use  of  (32)  in  (33)  enables  us  to  write 
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The  physical  interpretation  of  this  equation  is  as  follows.  The  spheroid  has  been  divided  up  to  M  belts, 
each  of  surface  A.v  =  7.7ili',li"An  ,  as  shown  in  Figure  2,  with  the  NSMS  being  an  unknown  constant 

over  each  belt.  At  the  center  of  each  segment,  the  sum  of  the  scattered  fields  from  all  M  belts  is  set  to 
equal  the  measured  field  Hdata(r)  at  point  r  that  is  a  known  field  arising  from  the  scatterer.  For  a  point  r„ 
the  latter  equation  leads  to 

M 

Tfl  f()|  ,r)  =  H>)  (35) 

m  v  'm  n 7  v  n7  v  7 

m= 1 

So  far  we  have  only  generated  one  equation  (or  three  if  we  have  access  to  the  full  vector  field)  with  M 
unknowns.  We  can  obtain  additional  independent  equations  by  using  data  collected  at  different  points  r„ 
with  n=  1,2,  N.  Matching  the  modeled  scattered  magnetic  field  to  the  data  at  these  N points  results  in 
the  linear  system 

[ZJ[QJ  =  [Hdata(r)L  (36) 

with 

"f(7,,r,)  f(/72,rj)  •••  f(77M,r,) 

\7  1=  f(?7''r2)  f^r2)  •••  f07 M’r2) 

_f(7,,rw)  f(tj2,rN)  ■■■ 

a  -nM]T,  (37) 

[Hdam(r„)]  =  [Hdata(r1)  Hda,a(r2)  Hdata(r;v)]T , 


and  r„)  given  by  (34),  whose  solution  can  be  written  symbolically  as 


[Q  ]  = 

L  in J 


[ZJT[Hdala(r  )] 
[ZJT[  ZJ 


(38) 


Once  [Q„,]  is  determined  the  object’s  EMI  response  can  be  computed  readily.  The  resulting  discrete 
NSMS  distribution  can  then  be  used  to  compute  the  total  NSMS  amplitude,  which  is  a  global  measure  of 
Q  for  the  entire  object  and  can  be  used  for  discrimination: 

M 

Q  =  Yn  AS  .  (39) 

^  mm  v  7 

m=l 

In  the  following  sections  we  study  some  features  of  this  global  measure  of  response. 
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2.3.4  The  dipole  model  as  a  limiting  case  of  NSMS 

Here  we  show  that  NSMS  reduces  in  the  limit  to  the  point  dipole  model  [6-11]  of  Section  2.2. 
Recall  that  the  magnetic  field  due  to  a  dipole  of  moment  m  is 

H“(r)  =  — Uj- m  (3RR-I),  (40) 

4;rR3 

where  R  is  the  unit  vector  along  R  =  r  -  rd  and  rd  and  r  are  respectively  the  location  of  the  dipole  and  the 

observation  point,  as  seen  in  Figure  3,  while  I  is  the  identity  dyad.  The  relation  between  the  induced 
dipole  moment  m  and  the  primary  magnetic  field  Hpr  at  the  dipole  location  is  given  by 

m  =  1^1  •  Hpr  (rrf ) ,  (41) 


where  the  magnetic  polarizability  tensor  M  depends  on  the  scatterer’s  shape,  size,  and  material 
properties.  For  a  body  of  revolution,  the  polarizability  tensor  in  a  coordinate  system  aligned  with  the 
scatterer’s  principal  axes  can  be  written  as 


M  = 


P  0 

'  no 


PP 

0  p  0 

r  pP 

oo/? 


(42) 


where  the  degeneracy  in  the  “radial”  element  Ppp  displays  the  BOR  symmetry  explicitly.  The  target’s 
principal  axes  and  the  global  coordinate  system  are  related  by  the  Euler  rotation  tensor. 

Now  let  us  prove  that  in  the  dipole  model  is  a  limited  case  of  the  NSMS.  To  do  that,  first  let  us 
divide  the  surrounding  spheroidal  surface  into  three  belts  and  assume  that  on  the  m-th  belt  the  NSMS 
density  follows  a  Dirac  delta  distribution  (see  Figure  3).  With  these  assumptions  the  scattered  magnetic 
field  (40)  becomes 

Hsc(r)  =  •  (3RmRm  - T)  (43) 

4  71  m=i  Rm 

where  now  R  =  r  -  rm  points  from  rm  on  the  m-th  belt  to  the  observation  point.  As  S  — >  0  we  have  that 
r,„  — >  rd,  and  ///'  ( rm )  =  Hf"'  ( r, ) ,  and  because  n2  =  p  =  cr,x  +  <x,y ,  and  n,=-n3=z,  ( a}  —  cos  a  and 

a,  =  sin  a  ,  where  a  is  the  angle  between  p  and  x )  then  (43)  reduces  to 

2DD  _  T 

Hsc(r)  =  — — •(2Q1//f(rd)z  +  Q2(//pr(rd)«1x  +  //f(rrf)a2y))  (44) 

4  nR 
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in  terms  of  the  Cartesian  unit  vectors  x ,  y  ,  and  z  .  After  introducing  a  diagonal  tensor 

~n2  o  o 

Mb=  0  n2  0  (45) 

0  0  2Q, 
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2.3.5  Invariance  properties  of  the  total  NSMS 

One  of  the  distinguishing  features  of  the  NSMS  model  is  that  the  total  NSMS  amplitude  is  a 
global  characteristic  of  the  object  and  is  independent  of  the  size  of  the  surrounding  spheroidal  surface.  To 
illustrate  this  feature  let  us  consider  a  sphere  of  radius  a,  conductivity  a,  and  permeability  //  =  jUrJuQ 
illuminated  with  a  time-varying  uniform  magnetic  field  Hpr  =  z  .  The  scattered  magnetic  field  at  any  point 
outside  the  sphere  can  be  written  in  analytic  form  as 

(47) 


where  ms  h  =  /7Hpr  is  the  magnetic  dipole  moment  of  the  sphere,  and  the  polarizability  /?  can  be  expressed 
in  terms  of  the  induction  number  k  =  ^jcou^up  through 


,  (2 //  + 1  )(ka coth ka-\)~ (kaf 

P  =  2 na  - - —  .  (48) 

( ju  - 1  )(ka  coth  ka-  \)  +  (ka) 


The  scattered  field  due  to  the  NSMS  density  in  the  spherical  coordinate  system  is 

Hsc(r)  =  V  [  ^ a1  dcpV D.{0) cos d^—^&inddd ,  (49) 

JO  JO 


where 


ii'  =  r'  =  sin0cos<px  +  sint?sin<py  +  cos  9z, 
h'  R  =  n' •  (r - r)  =  zQ cos 9 -a, 

R  =  ^ fa 2  +  Zq  -  2 az0  cos  9. 


(50) 


For  a  homogeneous  sphere  the  amplitude  of  NSMS  density  Q  =  Qq  is  constant  over  the  spherical  surface 
and  is  the  only  unknown  to  be  determined.  This  can  be  achieved  by  matching  the  secondary  magnetic 
field  of  the  NSMS  with  the  dipole  field  (40)  at  an  observation  point;  for  convenience  we  take  r  =  zQ z  . 
The  closed-form  solution  for  the  sphere  becomes 


H*ph  = 


m 

_ Z_ 

4/r 


2m 

_ Z 


\nz\ 


(51) 


and  the  scattered  NSMS  magnetic  field  due  to  a  uniform  source  density  Q(l  simplifies  to 
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(53) 


The  total  NSMS  is  thus  independent  of  the  size  of  the  surrounding  surface  and  is  an  intrinsic  property  of 
the  scatterer. 


2.3.6  Interpretation  of  the  total  NSMS 

The  total  NSMS  (and  its  time  evolution)  depends  on  the  size,  geometry,  and  material  composition 
of  the  object  in  question.  Early  time  gates  bring  out  the  high-frequency  response  to  the  shutdown  of  the 
exciting  field;  the  induced  eddy  currents  in  this  range  are  superficial,  and  a  large  NSMS  amplitude  at 
early  times  correlates  with  large  objects  whose  surface  stretches  wide.  At  late  times,  where  the  eddy 
currents  have  diffused  completely  into  the  object  and  low-frequency  harmonics  dominate,  the  EMI 
response  relates  to  the  metal  content  (i.e.,  the  volume)  of  the  target.  Thus  a  smaller  but  compact  object  has 
a  relatively  weak  early  response  that  dies  down  slowly,  while  a  large  but  thin  or  hollow  object  has  a 
strong  initial  response  that  decays  quickly.  These  features  can  be  neatly  summarized  by  the  parameters  of 
an  empirical  decay-law  model  like  the  Pasion-Oldenburg  law  see  (57). 

2.3.7  The  parameterized  NSMS 

During  APG  standardized  test-site  discrimination  studies  (see  Chapter  6)  we  use  a  parameterized 
version  of  NSMS  to  encapsulate  the  electromagnetic  signature  of  a  target  [14].  In  this  version  of  the 
model — which  provides  at  least  three  independent  polarizability-like  parameters  for  use  in  discrimination 
and  thus  in  a  sense  extracts  further  information  from  the  same  data — the  scatterer  is  associated  with  a 
surrounding  sphere  S  on  which  a  set  of  dipoles  are  distributed.  The  secondary  field  is  expressed  as 

~  ~  |Tn„(r,„0  r^r(«»' 

Hsc(r,0  =  4  - n„(r,,0  nw(r ,,t)  nyz(rs„t )  H* Ov)  ds'  =  Z-Cl,  (54) 

4-7rKc, 

[ac(rs„0  nzy(r  ,,t)  T>  „ (rs, , 0 JJ  \_Hzr (rs, ) 
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where  R  ,  points  from  the  location  r  of  the  s'  -th  patch  on  the  sphere  to  the  observation  point  r  and  the 
response  amplitude  of  each  patch  is  a  combination  of  the  primary  field  piercing  it  and  the  tensor  of 
normalized  strengths  Q  (r,,t),  which,  as  usual  [15],  is  symmetric:  Q  =Q  .  The  z  -axis  is  dictated  by 
the  direction  of  m  from  HAP  or  from  the  dipole  model,  and  the  x  -  and  y  -axes  are  arbitrarily  chosen  to 

be  perpendicular  to  z  and  to  each  other.  The  integral  is  again  transformed  to  a  matrix-vector  product 
through  numerical  quadrature.  The  amplitude  array  Q  is  determined  by  minimizing  in  a  least-squares 
sense  the  difference  between  measured  data  with  a  known  object-sensor  configuration  and  the  predictions 

of  equation  (54).  Once  the  tensor  elements  Q.  (s')  are  found  one  can  define  “total  polarizabilities”  by 

y 

integrating  over  the  sphere, 

Qij(t)  =  lni.(r„t)ds',  (55) 

and  these  can  in  turn  be  used  to  find  “principal  elements”  through  joint  diagonalization: 


QJ  0 

QJ  0 

2,(0 

"  2,(0 

0 

0 

2,(0 

<2,(0 

QJ  o 

=  A 

0 

2V(0 

0 

QJ  0 

2,(0 

2,(0 

0 

0 

2,(0  _ 

where  the  matrix  A  is  orthogonal  and  the  prime  denotes  transposition.  The  information  contained  in  the 
diagonal  tensor  can  be  summarized  further  by  incorporating  the  empirical  decay  law  of  Pasion  and 
Oldenburg  [16]: 

MJt)  =  Qcp)  =  Bjp««e7°J,  a  =  x,  y,  z,  (57) 

where  t  is  the  time,  B  ,  j3aa  ,  and  yaa  are  the  fitting  parameters,  and  Maa(t )  is  the  total  NSMS  along 

the  x,  y,  and  z  directions  in  the  body  frame.  The  principal  NSMS  elements  and  the  Pasion-Oldenburg 
parameters  are  intrinsic  to  the  object  and  can  be  used,  on  their  own  or  in  combination  with  other 
quantities,  in  discrimination  processing. 

2.4  The  orthonormalized  volume  magnetic  source  model 

Most  EMI  sensors  are  composed  of  separate  transmitting  and  receiving  coils.  When  the  operator 
activates  the  sensor,  a  current  runs  through  the  transmitter  coils,  which  results  in  the  establishment  of  a 
(“primary”  or  “principal”)  magnetic  field  in  the  surrounding  space  (Figure  4).  According  to  the 
elementary  atomic  model  of  matter,  all  materials  are  composed  of  atoms,  each  with  a  positively  charged 
nucleus  and  a  number  of  orbiting  negatively  charged  electrons.  The  orbiting  electrons  cause  circulating 
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currents  and  form  microscopic  magnetic  dipoles.  In  the  absence  of  an  external  magnetic  field  the 
magnetic  dipoles  of  atoms  of  most  materials  have  random  orientations,  resulting  in  no  magnetic  moment. 
The  application  of  an  external  time  varying  magnetic  field,  by  Faraday’s  law,  induces  eddy  currents  in 
highly  conducting  bodies  by  an  alignment  of  the  magnetic  moments  of  the  spinning  electrons  and  a 
magnetic  moment  due  to  a  change  in  the  orbital  motion  of  electrons.  These  currents  and  magnetization  in 
turn  generate  a  (“secondary”  or  “scattered”)  magnetic  field  that  also  varies  with  time  and  induces 
measurable  currents  in  the  receiving  coils.  The  induced  magnetic  dipoles/eddy  currents  are  distributed 
inside  the  object  and  produce  a  magnetic  field  intensity  H  outside.  The  magnetic  field  due  to  the  i  -th 
source  can  then  be  expressed  at  any  observation  point  r  as  the  matrix-vector  product 

H.(r)  =  G.(r)m.,  (58) 


where  the  Green  function  G.  is  given  in  detail  in  equation  (1).  When  there  are  several  such  sources,  the 
total  field  can  be  expressed  as  a  superposition: 


H(r)  =  £G,(r)m,= 


G,  G2 


m 


m 


(59) 


Before  going  further  we  note  that  our  method  takes  as  input  the  (in  principle  unknown)  number 
M  of  radiating  sources.  For  advanced  EMI  sensors  such  as  the  MetalMapper  and  2x2  and  5x5 
TEMTADS  arrays  we  have  developed  a  procedure  based  on  joint  diagonalization,  sketched  in 
Section  2.5,  that  estimates  M  stalling  from  raw  data  and  with  no  need  for  inversion.  For  other  sensors 
one  may  proceed  by  letting  M  vary  as  part  of  an  optimization  routine. 


Figure  4:  A  metallic  object  under  the  transmitter.  The  target’s  EMI  response  at  the  receiver  coil  can  be  calculated 
from  the  equivalent  surface  or  volume  magnetic  dipole  moment  dm. 
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The  superposition  (59)  can  be  used  (and  often  has)  to  carry  out  one-  and  multi-object  inversions 
starting  from  data  taken  at  an  ensemble  of  points.  All  the  measured  H  -values — which  can  pertain  to 
multiple  transmitters,  multiple  receivers,  and  different  vector  components — are  strung  together  in  a  one¬ 
dimensional  array,  while  the  corresponding  Green  functions  are  stacked  as  matrix  rows.  The  resulting 
composite  G  matrix  can  then  be  (pseudo)inverted  to  find  the  strengths  of  the  sources.  This  procedure, 
which  is  nothing  other  than  the  dipole  model  if  each  body  is  taken  to  be  represented  by  one  source  only, 
works  well  for  one  or  two  sources,  but  for  larger  numbers  becomes  very  time-consuming  (since  the  Green 
matrix  becomes  very  large)  and  increasingly  ill-posed,  usually  requiring  regularization.  The  ONVMS 
method  is  designed  to  circumvent  these  difficulties. 

2.4.1  Orthonormal  Green  functions 

The  method  starts  from  the  realization  that  the  matrix-vector  product  (58)  is  valid  at  any 
observation  point  r  and,  in  particular,  at  every  point  r  .  If  we  introduce  the  inner  product 

(a.b)=[  ATBcls=[  ATBds+\  ATBds  +  ...,  (60) 

'  '  J  S  j  Rx0  J  Rx, 

where  the  integral  is  computed  over  the  “sensitive”  surfaces  of  the  sensor,  and  if  furthermore  we  can  find 
a  basis  of  Green  functions  orthogonal  under  this  measure, 

M 

H(r,)  =  2>/r)b  such  that  (¥,,*,)  =  F/,,,  (61) 

7=1 

where  S jk  is  a  Kronecker  delta,  then  it  is  possible  to  find  the  source  amplitudes  b  .  without  costly  and  ill- 
conditioned  inversions  simply  by  exploiting  the  sifting  property  of  the  orthogonal  basis: 

M  M 

7=1  7=1 

and  thus 

(63) 

which  clearly  does  not  involve  solving  a  linear  system  of  equations;  it  is  necessary  to  invert  only  the  6x6 
matrix  F  .  Moreover,  this  definition  of  the  coefficients  b  guarantees  that  they  are  “optimal”  in  the  sense 

k  J 

that  the  expansion  (61)  yields  the  least  mean-square  error  ^H-  X^TVb^H-  X^TVb.^  [17]. 
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To  construct  the  set  of  orthonormal  Green  functions  we  resort  to  a  generalization  of  the  Gram- 
Schmidt  procedure  [18].  Assuming  that  the  Green  matrices  are  linearly  independent — i.e.,  that  we  cannot 
have  a  collection  of  distinctly  located  dipole  sources  combining  to  produce  no  measurable  field  unless 
their  amplitudes  all  vanish — we  define 

*i  =  <V 

*2  =  G2-¥Ai. 


m—  1 

¥  =G  -V TEA  (64) 

m  m  k  mk  7 

k= 1 

M- 1 

—  n  —  V  m  a 

M  L-i  k  Mk’ 

k= 1 

where  the  6x6  matrices  A  obey  A  - 0  for j<k.  Enforcing  the  orthogonality  relation  (61)  is 
equivalent  to  setting  (^¥n,Gm}  =  FnAm  for  n  <  rn ,  and  using  this  relation  twice  in  definition  (64)  we  find 

f  "-1  ^ 

Ann  =  FA  ^  C n,n  ~  Z  Kk  A  A mk  J>  (65) 

where  the  overlap  integral  C  =(G  ,G  ). 

A  mn  \  m  n  / 

At  the  end  of  the  process  it  is  necessary  to  recover  an  expansion  expressed,  like  (58),  in  terms  of 
the  actual  Green  functions,  in  part  because  the  functions  VF  .  are  orthogonal  (and  defined)  only  at  points 

on  the  receivers,  and  in  part  because  of  the  non-uniqueness  of  the  coefficients  b  due  to  the  arbitrary 
order  in  which  the  G  enter  the  recursion  (64).  To  that  end,  we  express 


m  k  m 


(66) 


and  to  find  the  coefficients  B  we  compare  expansion  (66)  term  by  term  to  the  definition  (64)  and  use 
the  rule  that  A  =  0  for  j  <k  to  find 


B  -I,  the  identity, 

mm  7  J  7 


R  =—A 

m(m-l)  m(m-l)  ’ 


B  -  -V  B.  A  ,  for  1  <  q  <  m  -  2, 

mq  Iq  ml  1 


i=q 


(67) 
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in  terms  of  which  we  recover  the  physical  polarizability  elements: 

M  M  f  k  A  M  f  M  Am 

H  =  2Xb,  =  Z|  I G,fl„  I  b,  =  Zg,  I  £uHb,  I  =  •  «»> 

k= 1  k= 1  V  /= 1  J  1= 1  V  k=l  J  1=1 

2.4.2  ONVMS  procedure 

With  all  the  pieces  in  place,  we  can  sketch  an  algorithm  to  invert  EMI  data  using  the  ONVMS 

model: 

1)  Given  a  number  of  sources  and  their  tentative  locations,  find  the  Green  tensors  G,  =T  using 
equation  (5)  and  compute  the  overlap  integrals  G  using  the  inner  product  (60). 

2)  Determine  the  first  normalization  factor,  F{  =  (CJ^G^ ,  and  use  it  to  find  all  the  Gram-Schmidt 
coefficients  A  with  n  =  1 :  A  =  F  C,  . 

mn  ml  1  1  m 

3)  Set  m  =  2 ;  compute,  in  sequence, 

a)  The  coefficients  Amn  with  n  =  2,  ,m-l  using  equation  (65); 

b)  The  function  xPj  using  the  expansion  (64); 

c)  The  normalization  factor  F  =/xF  ,  T  )  ; 

y  m  \  m  m  / 

increase  m  by  1  and  iterate  until  all  sources  have  been  included. 

4)  Once  all  the  A  ,  F  ,  and  T  are  known,  find  B  using  (67). 

5)  Use  the  orthonormality  of  the  new  Green  functions  to  determine  the  source  amplitudes  using 
b  =  G  1  ^xF?,Hdata^  ,  as  in  (63).  Take  the  measured  field  to  be  piecewise  constant — i.e.,  constant 
throughout  each  receiver — when  evaluating  the  integrals. 

6)  Use  the  computed  ,  B  ,  and  Gm,  along  with  the  expansion  (68),  to  generate  the  secondary 
field  prescribed  by  the  given  number  of  sources  at  the  given  locations. 

7)  Compare  the  model  prediction  with  the  measured  data,  vary  the  source  locations,  and  iterate  until 
the  least-squares  discrepancy  between  prediction  and  measurement  attains  a  suitable  minimum. 

The  procedure  as  written  applies  to  only  one  time  gate,  but  the  extension  to  fully  time-dependent 
functions  is  straightforward:  we  need  only  substitute  the  vectors  b  and  Hdata  for  two-dimensional  arrays 
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where  the  columns  denote  time.  The  relations  between  the  two,  namely  (63)  and  (68),  acquire  multiple 
right-hand-sides,  and  the  optimization  mentioned  on  Step  7  of  the  algorithm  is  constrained  further.  As  a 
final  remark  we  note  that  rigorously  speaking  the  coefficients  b  (and,  for  that  matter,  the  amplitudes 

q 

m  )  are  not  the  polarizabilities  themselves  but  relate  more  closely  to  their  time  derivatives  [19-20]. 

The  great  advantage  of  the  ONVMS  technique  is  that  it  takes  into  account  mutual  couplings 
between  different  parts  of  targets  and  avoids  matrix  singularity  problems  in  cases  with  multiple  objects. 
Once  the  polarizability  tensor  elements  and  the  locations  of  the  elemental  responding  dipoles  are 
determined  one  can  group  them  according  to  their  volume  distribution.  For  each  group  a  total 
polarizability  tensor  can  be  computed  and  diagonalized  using  joint  diagonalization,  the  topic  of 
Section  2.5.  The  resulting  time -dependent  diagonal  elements  have  been  shown  to  be  intrinsic  to  the 
objects  and  can  be  used,  on  their  own  or  combined  with  other  quantities,  in  discrimination  processing. 

2.5  Joint  diagonalization  for  multi-target  data  pre-processing 

In  real  life  situations  the  targets  of  interest  are  usually  surrounded  by  natural  and  artificial  debris 
with  metallic  content,  including,  for  instance,  the  remains  of  ordnance  that  did  explode.  Thus  it  is  usually 
not  clear  how  many  objects  are  producing  a  given  detected  signal;  all  sensing  methods,  including  EMI, 
are  fraught  with  detection  rates  that  overwhelm  cleanup  efforts  and  hike  their  cost.  Here  we  introduce  a 
data  pre-processing  technique  based  on  joint  diagonalization  (JD)  that  estimates  the  number  of  targets 
present  in  the  field  of  view  of  the  sensor  as  it  takes  a  data  shot,  and,  in  a  good  number  of  cases,  even 
provides  the  capability  to  perform  real-time  characterization  and  classification  of  the  targets  without  the 
need  for  a  forward  model. 

Joint  diagonalization  has  become  an  important  tool  for  signal  processing  and  inverse  problems, 
used  as  part  of  independent  component  analysis  [21],  blind  source  separation  or  BSS  [22],  common 
principal  component  analysis,  and,  more  recently,  kernel-based  nonlinear  BSS  [23].  We  further  extend  the 
applicability  of  the  method  by  using  it  to  detect  and  locate  buried  targets  without  the  need  for  inversion. 
As  we  say  above,  a  variation  of  the  method  can  be  used  to  extricate  time -dependent  electromagnetic 
signatures  from  attitude  information.  Here  we  will  outline  the  detailed  procedure  as  applied  to  the 
TEMTADS  sensor  array,  a  time-domain  device  with  25  transmitter/receiver  pairs  that  provides  625 
measurements  over  Ng  =  123  time  gates  at  each  sensor  location. 

2.5.1  The  multi-static  response  matrix 

JD  estimates  the  eigenvalues  and  eigenvectors  of  a  square  time-  or  frequency-dependent  multi¬ 
static  response  (MSR)  matrix  synthesized  directly  from  measured  values.  To  construct  the  MSR  matrices 
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one  just  has  to  stack  the  625  readings  at  each  time  gate  in  a  25  x  25  array  so  that  each  column  stands  for 
one  of  N,  transmitters  and  each  row  represents  one  of  Nr  receivers: 


,  k  =  l,K,Ng,  (69) 

where  the  element  Hu  is  the  field  measured  by  the  i-th  receiver  when  the  j-th  transmitter  is  fired.  The 
second  step  of  the  procedure  is  to  diagonalize  the  123  matrices  at  one  stroke  so  they  all  share  a  single  set 
of  orthonormal  eigenvectors.  In  other  words,  given  the  MSR  matrix  S (tk)  at  the  k- th  time  gate,  we  look  for 
a  unitary  matrix  V  such  that  the  products 

D*=VrS(gV  (70) 

are  “as  diagonal  as  possible”  (i.e.,  their  off-diagonal  elements  vanish  within  a  preset  tolerance).  By 
diagonalizing  all  the  matrices  simultaneously  we  separate  the  time -dependent  intrinsic  features  of  the 
responding  sources  (and  hence  the  interred  objects),  which  get  encapsulated  in  the  eigenvalues,  from  the 
other  factors — notably  the  location  and  orientation  of  the  target  with  respect  to  the  sensor — that  influence 
the  signal  but  do  not  change  as  the  data  are  being  taken;  these  get  bundled  into  the  eigenvectors.  (The  fact 
that  the  locations  and  orientations  can  be  dissociated  in  this  way  from  the  electromagnetic  signatures  is  an 
upside  of  the  low  frequencies  of  the  quasistatic  EMI  range,  because  the  relevant  Green  functions  are  time- 
independent.)  Thus  the  measured  data  can  be  resolved  as  a  supeiposition  of  “elemental”  sub-signals,  each 
corresponding  to  an  elementary  dipolar  source,  whose  combination  corresponds  to  the  buried  objects. 
Each  source — and  the  corresponding  field  singularity — can  moreover  be  localized  numerically:  the 
TEMTADS  geometry  is  such  that  the  diagonal  of  the  unprocessed  MSR  matrix  mimics  a  set  of 
monostatic  measurements,  akin  to  those  taken  with  a  handheld  sensor,  which  peak  sharply  when  there  is  a 
target  directly  underneath.  The  maxima  in  the  diagonal  thus  point  to  the  transmitter/receiver  pairs  closest 
to  any  responding  sources.  These  location  estimates  can  be  grouped  and  correlated  to  the  eigenvalue 
distributions  to  estimate  target  locations. 


H„ 

7/ 

11 

12 

IN, 

//  , 

// 

/A 

21 

22 

2  N, 

M 

M 

0 

M 

HN2 

L 

Hnn, 

2.5.2  Interpretation  and  diagonalization  of  the  MSR  matrix 

We  now  proceed  to  express  our  above  considerations  quantitatively.  Initially  we  consider  the 
transmitter  assembly,  which  in  TEMTADS  consists  of  a  set  of  coplanar  square  loops  forming  a  regular 
grid.  The  Biot-Savart  law  gives  the  primary  magnetic  induction  established  at  the  location  r,  of  the  /-th 
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source  when  the  j- th  transmitter  antenna  (whose  area  is  <rTx  )  is  excited  immediately  before  shutoff  by  a 
current  If 


Bp,r  = 

ji 


Mo1! 
4  n 


dl'  x  (r  -r') 


lr,-r 


'  |3 


=  ^Jr 


(71) 


This  primary  field  induces  in  the  /-th  source  a  dipole  moment  given  by 

m,  =  U,A,UfB£, 


(72) 


where  the  Euler  rotation  matrix  U  relates  the  instrument’s  coordinate  axes  to  the  principal  axes  of  the 
source,  and  the  diagonal  polarizability  matrix  A the  only  quantity  intrinsic  to  the  source,  measures  the 
strength  with  which  the  primary  field  induces  a  moment  along  each  of  those  axes. 

According  to  Faraday’s  law,  the  signal  measured  by  a  receiver  coil  is  the  electromotive  force 
given  by  the  negative  of  the  time  derivative  of  the  secondary  magnetic  flux  through  the  coil.  Since  the 
field  at  point  r  of  a  dipole  of  moment  m  placed  at  r0  is  given  by 

ti  c  j. _ j.  i  ti  j. _ j. 

B  =  —  Vxlmx - I ,  and  thus  [  B  •  ds  -  -m  •  —  Xafl  x - \  (73) 

4/r  ^  lr-r0l3J  J  4/r®  lr-r0l3 


by  straightforward  application  of  Stokes’s  theorem,  one  obtains  that  the  signal  sampled  at  time  tk  by  the  i- 
th  receiver  (of  area  crRx  )  when  the  /-th  source  is  excited  by  the  j- th  transmitter  is 


Mo 


Htj(.tk)<y rx,^tx/;  =  ^"rx 

1  '  4k 


1  ^'X(r'-r<> 


™j,(tk)  =  -m77 (fk) 

=  g,X,'[UA((4)U']-gX/p 


(74) 


where  a  dot  over  a  variable  indicates  its  time  derivative.  In  equations  (71)  and  (74)  the  line  element  dl' 
lies  on  the  x-y  plane,  and  as  a  consequence  the  Green  functions  are  similar  in  structure  to  those  of  the 
simple  model  presented  in  Section  2.2.  Note  that  we  have  included  the  exciting  current  /,  and  the 
transmitter  and  receiver  areas  in  the  definition  of  the  signal;  we  have  explicit  knowledge  of  these 
quantities  and  can  factor  them  out.  If  only  the  /-th  source  is  illuminated,  we  construct  the  MSR  matrix  for 
the  complete  transmitter/receiver  array  by  tiling  Nr  x  N,  instances  of  the  expression  (74): 

S  =  GscU,A/U^(Gpr)r,  (75) 
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where  the  primary  (or  transmitter)  dyad  Gpr  is  of  size  N,  x  3,  the  secondary  (or  receiver)  dyad  G  c  is  of 
size  Nr  x  3,  and  the  response  matrix  IJA/IJ7  is  3  x  3.  When  there  is  more  than  one  source  present,  the  MSR 
matrix  of  equation  (75)  is  readily  generalized: 


s  =  [g*c  G*c  ...] 


=  Gru,  Gru2  ...] 


"UjAjUf  0 

"(erf 

0  U2A2U2  - 

(GHr 
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>■ 
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_ * _ i 

'(Gru/' 

0  a2  ••• 

(GprU2)r 

_ 

(76) 


where  we  see  that  the  features  intrinsic  to  the  targets  can  be  separated  formally  from  the  particulars  of  the 
measurement — that  is,  from  the  geometry  and  dimensions  of  the  sensor  and  the  sensor-target  attitude.  The 
array  S  has  size  Nr  x  N,  and  is  square  if  Nr  =  N„  as  is  the  case  with  TEMTADS.  This  allows  us  to 
diagonalize  the  matrix  but  does  not  suffice  to  guarantee  that  the  extracted  information  is  useful — i.e.,  that 
the  eigenvalues  and  eigenvectors  are  real,  and  that  the  latter  are  orthonormal.  For  that  to  hold  we  must 

have  a  real,  symmetric  matrix,  which  requires  G“  =  Gpr  =  G; .  This  cannot  be  rigorously  true,  because 

the  receivers  cannot  coincide  exactly  with  the  transmitters,  but  holds  approximately  for  TEMTADS  if  we 
factor  the  exciting  current  and  the  coil  areas  out  of  S,  as  we  did  in  equation  (74).  The  diagonalization  we 
perform  is  thus  a  particular  case  of  a  singular  value  decomposition  (SVD),  and  in  what  follows  we  use 
“diagonalization”  as  shorthand  for  “SVD  of  a  symmetric  matrix.” 

The  decomposition  (76)  exhibits  the  actual  polarizability  elements  but  is  not  directly  available  to 
us  because  the  Green  tensors  are  not  orthogonal.  To  see  what  we  do  get  when  we  diagonalize  S  we  can 
perform  the  SVD  on  G: 

S  =  GUAUrGr  =  w[zVrUAUr  Vz]  Wr  =  WZAZrWr  =  YAYr  (77) 

In  the  intermediate  step  we  have  used  the  fact  that  the  matrix  within  the  brackets  is  real  and  symmetric 
and  thus  has  a  purely  real  eigendecomposition.  Result  (77)  shows  that  the  eigenvalue  matrix  A,  though 
time-dependent,  is  not  solely  composed  of  source  responses,  but  also  contains  location  and  orientation 
information  extracted  from  the  Green  tensors.  The  eigenvectors,  likewise,  include  information  from  both 
the  polarizabilities  and  the  measurement  particulars. 

We  also  see  in  the  decomposition  (77)  that  S  contains  an  unknown  “hidden  dimension” — 3 N, 
where  N  is  the  number  of  sources — in  the  size  of  the  block-diagonal  response  matrix.  Numerical 
diagonalization  (or,  in  general,  the  SVD)  of  S  will  impose  this  middle  dimension  to  be  Nr  =  N,.  Ideally, 
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the  method  should  be  able  to  resolve  up  to  ^Nr  /  3j  responding  sources,  or  eight  for  TEMTADS,  but  the 

actual  number  is  lower.  For  one,  the  procedure  will  resolve  targets  only  when  they  are  spatially  separated: 
two  distinct  dipoles  sharing  one  location  decrease  the  rank  of  the  G  matrices,  and  hence  of  S,  by  3.  In  any 
case,  diagonalization  of  S  can  again  let  us  estimate  the  number  of  targets  illuminated  by  the  sensor;  since 
the  only  time -dependent  quantities  are  the  intrinsic  polarizabilities  of  the  sources,  we  expect  the 
additional  information  provided  by  the  time  decay  of  the  eigenvalues  to  be  useful  for  classification. 

The  development  outlined  above  corresponds  to  each  time  gate  taken  separately.  To  make  sense 
of  the  time-dependent  information  we  have  to  find  a  way  to  “follow”  each  of  the  eigenvalues  as  the  signal 
decays.  (A  similar  process  must  be  carried  out  when  using  the  dipole  model  for  inversion.)  One  could  in 
principle  diagonalize  the  MSR  matrix  at  each  time  channel,  and  the  eigenvectors,  which  depend  only  on 
geometry  and  pose,  should  stay  constant;  however,  it  is  not  possible  to  know  a  priori  the  order  in  which 
the  eigenvalues  will  be  given  by  the  diagonalization;  this  fact — not  to  mention  noise  and  experimental 
uncertainty — makes  it  inevitable  to  have  to  disentangle  the  tensor  elements  by  hand,  which  is  easily  done 
wrong.  Instead,  we  explicitly  look  for  an  orthogonal  matrix  of  eigenvectors  that  diagonalizes  all  the  MSR 
matrices  simultaneously.  The  procedure  we  employ  is  a  generalization  of  the  method  for  single  matrices, 
and  is  well-known;  it  is  sketched  in  next  Section. 

2.5.3  Algorithm  for  joint  diagonalization 

The  joint  diagonalization  algorithm  we  use  [22,  24-25]  is  a  generalization  of  Jacobi’s  procedure 
to  find  the  eigenvalues  of  a  single  matrix.  Formally  we  set  out  to  solve  the  optimization  problem 

(78) 

avA(tqwTij)2 

^  q= 1  i*j 

s.t.  VTV  =  /, 

which  we  accomplish  by  making  repeated  Givens-Jacobi  similarity  transformations  designed  to  gradually 
accumulate  the  “content”  of  the  matrices  on  their  diagonals  until  a  certain  tolerance  level  is  reached.  The 
transformations  are  of  the  form  A(t  )  — >  A'(t  )  —  VmA(t  )V* ,  with  the  matrix  Vrs  being  the  identity  but 
with  the  four  elements  Vrr,  Vrs,  Vsr,  and  Vss  replaced  by  the  two-dimensional  rotation  array 

with  tan  2</>rs  = - frs  ,  (79) 

n  +  Jf2  +  n 2 

rs  V  J  rs  fs 

where 


cos  (b  sin  d) 

r  rs  '  rs 

-  sin  (b  cos  d) 
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n„  =X{K-(*»)-fl®(*»)]2 -Mtq)  +  dsr(tq)f},  (80) 

frs  =2^[arr(tq)-aa(tq)][als(tq)  +  asr(tq)].  (81) 

The  indices  are  swept  systematically,  and  the  procedure  is  repeated  until  convergence  is  reached. 
The  computational  burden  is  equivalent  to  that  of  diagonalizing  the  matrices  one  by  one.  The  resulting 
eigenvalues  and  eigenvectors  are  all  real  because  all  the  MSR  matrices  are  symmetric. 
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3  Inverse  Models 

3.1  Introduction 

Several  EMI  sensing  and  data-processing  techniques  [1-4,  7,  13,  19,  26-38]  have  been  recently 
developed  for  detecting  and  discriminating  between  UXO  and  non-UXO  items.  Typically  the  first  step  of 
these  methods  is  the  recovery  of  a  set  of  parameters  that  specify  a  physics-based  model  representing  the 
object  under  interrogation.  For  example,  in  EMI  sensing,  the  recovered  parameters  consist  of  the  object’s 
location  and  spatial  orientation  in  addition  to  “intrinsic”  parameters  such  as  the  polarizability  tensor 
(along  with  some  parameterization  of  its  time -decay  curve)  in  dipole  models  or  the  amplitudes  of 
responding  magnetic  sources  in  the  NSMS  and  ONVMS  models.  EMI  responses  depend  nonlinearly  on 
the  subsurface  object’s  location  and  orientation,  therefore  determining  the  buried  object’s  orientation  and 
location  is  a  non-linear  problem.  In  this  section  several  inverse  scattering  approaches  are  described  for 
EMI  data  inversion. 

Most  EMI  sensors  are  composed  of  separate  transmitting  and  receiving  coils.  When  the  operator 
activates  the  sensor,  a  current  runs  through  the  transmitter  coils,  resulting  in  the  establishment  of  a 
(“primary”  or  “principal”)  magnetic  field  in  the  surrounding  space.  By  Faraday’s  law,  this  time-varying 
magnetic  field  induces  eddy  currents  in  highly  conducting  bodies  (ferromagnetic  bodies  also  have  their 
magnetization  affected  by  the  impinging  field).  These  currents  and  magnetization  in  turn  generate  a 
(“secondary”  or  “scattered”)  magnetic  field  that  also  varies  with  time  and  induces  measurable  currents  in 
the  receiving  coils.  At  the  end,  the  electromagnetic  data  are  inverted  using  different  forward  models.  The 
procedure  for  estimating  the  location,  orientation,  and  electromagnetic  parameters  of  a  buried  object 
(linked  in  a  “model  vector”  v)  is  carried  out  by  defining  an  objective  function  that  quantifies  the 
goodness-of-fit  between  the  measured  data  and  the  predictions  of  the  forward  model.  Routinely,  a  least- 

squares  (LS)  approach  is  taken  to  recover  v:  formally,  if  dobs  is  the  vector  of  the  measured  scattered  field 
and  F(v)  is  the  solution  to  the  forward  problem,  the  least-squares  criterion  assumes  the  form 

minimize  cj){  v)  =  II  dobs  -  F(v)  II  .  (82) 


A  simple  way  to  determine  the  model  vector  v  is  to  use  the  Gauss-Newton  method,  which  starts 
with  an  initial  guess  v0  and  updates  it  iteratively  through 


v,  ,  =  v ,■ 

£+1  k 


(83) 


where  k  denotes  the  iteration  number  and  s *  is  a  perturbation  direction;  we  solve  for  the  S/  that  minimizes 
(/).  In  many  cases  the  LS  approaches  suffer  from  an  abundance  of  local  minima  that  often  leads  them  to 
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make  incorrect  predictions  of  location  and  orientation.  Global  search  procedures,  such  as  differential 
evolution  (DE)  [33-34]  and  genetic  algorithms  [37],  have  been  recently  developed  to  avoid  this  problem. 
We  have  combined  the  DE  algorithm  with  the  NSMS  model  [2]  (or  with  the  dipole  model  [37])  to  recover 
locations  and  orientations  of  buried  objects.  Once  these  extrinsic  properties  are  found  we  perform 
classification  using  Mixed  Models  (MM)  and  standard  Matlab  built-in  classifiers  based  on  maximum 
likelihood  methods  or  on  linear,  quadratic,  or  Mahalanobis  distances.  Both  gradient  and  global  search 
approaches  are  computationally  intensive  because  they  require  a  massive  number  of  forward-model 
evaluations  and  because  the  determination  of  the  nonlinear  elements  of  v — the  location  and  orientation  of 
the  object — is  a  nontrivial  and  time-consuming  problem  in  itself.  To  avoid  non-linear,  time-consuming 
inversions,  and  by  so  doing  streamline  the  inversion  process,  we  recently  developed  a  new  physics-based 
approach  called  the  HAP  method  and  applied  it  to  various  UXO  discrimination  problems.  The  HAP 
method  exploits  an  analytic  relationship  between  the  magnetic  field  vector  H ,  the  vector  potential  A , 
and  the  scalar  magnetic  potential  y/  (Psi)  of  a  hypothetical  point  dipole  to  determine  the  location  of  a 
visually  obscured  object.  Of  these  quantities  only  the  magnetic  field  (and  often  only  one  of  its 
components)  is  available,  and  as  part  of  this  project  we  developed  a  numerical  procedure  based  on  the  2D 
NSMS  model  that  replaces  the  measurement  surface  around  the  scatterer  with  a  flat  plane  of  dipoles  at  a 
(known)  location  intermediate  between  the  instrument  and  the  target.  The  amplitudes  of  these  responding 
sources  can  be  computed  starting  from  high-spatial-coverage  geophysical  data  by  solving  a  linear  system 
of  equations  and  can  then  be  used  to  reconstruct  H  ,  A ,  and  y/  at  any  point  on  or  above  the  measurement 

surface  and  thus  to  solve  for  the  relative  location  R  and  the  polarizability  M  of  the  hypothetical  dipole. 

This  chapter  briefly  overviews  gradient-based  optimization,  differential  evolution,  and  the  HAP 

method. 

3.2  Gradient-based  methods  of  optimization 

One  of  the  most  popular  approaches  for  solving  inverse  problems  is  the  gradient  method  [39-41], 
The  gradient  method  requires  the  system’s  Jacobian,  which  contains  the  gradients  of  the  scattered  field 
with  respect  to  the  unknown  parameters  of  interest.  In  many  cases  it  is  impossible  to  determine  the 
scattered  EM  field’s  derivatives  analytically;  this,  however,  is  not  a  problem  with  either  the  dipole  model 
or  the  NSMS  model.  Further,  the  NSMS-based  inverse  approach  always  results  in  an  over-determined 
system  and  thus  does  not  suffer  from  the  ill-conditioning  that  usually  afflicts  finite -element  or  finite- 
difference  time -domain  methods.  The  EM  scattering  problem  can  be  written  in  compact  matrix  form  as: 

[Z][Q]  =  {//"}  (84) 
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where  [Z]  is  the  scattering  matrix,  {Q}  is  a  vector  containing  the  amplitudes  of  responding  dipoles 
(normalized  by  the  primary  field),  and  {Hd }  is  a  vector  containing  the  measured  data  over  a  set  of  points. 
The  important  point  to  note  is  that  [Z]  in  the  NSMS  contains  explicit  expressions  for  the  responding 
source  amplitudes  {Q}  in  terms  of  the  object’s  location  and  orientation  that  can  be  differentiated 
analytically  and  that  contain  no  singularities  in  the  regions  where  they  must  be  evaluated.  Let  us  assume 
that  a  is  a  set  of  parameters  (orientation,  depth,  etc.)  that  must  be  determined  from  a  set  of  measured  data 
[42] .  A  convenient  way  to  view  the  problem  is  to  define  a  forward  map  as  one  that  associates  a  given  a 
with  an  initial  value  a0  (which  serves  to  kick-start  the  inversion  process).  A  least-squares  formulation  of 
the  problem  identifies  a  minimum  of  the  error  function  by  solution  of  the  equation 

KM7!,.,  <85> 

where  [/]a_,  is  a  Jacobian  matrix  based  on  { a/(  ( } .  ft  is  the  iteration  number,  the  modeled  values 
{ H "i°d { cc ]  ]  are  predicted  based  on  {a  },  and  the  solution  {Sap}  is  a  vector  of  incremental  steps  in 
the  unknown  parameters,  which  are  updated  via 

{afi}  =  {ap_1}  +  {Safi\.  (86) 

For  a  case  in  which  the  EMI  response  from  a  body  of  revolution  (BOR)  is  approximated  using 
five  NSMS  sources,  the  solution  vector  contains  10  numbers:  five  describe  its  location  and  orientation, 
0  =  (x  ,  y  ,  z0 ,  0,  cp) ,  and  five  more  are  the  source  amplitudes  (“omega  parameters”).  These  parameters  are 

recovered  from  measurements  H±ia  of  the  secondary  magnetic  field  at  the  center  of  the  sensor. 
Specifically,  the  parameters  are  found  from  minimization  of  the  sum  of  squares  (SS), 

1  Axp 

SS  =  -^(//. -//data)2  (87) 

2  (=i 

where  H. ,  the  theoretical  scattered  magnetic  field  at  each  of  JVexp  measurement  points,  is  a  function  of  10 
parameters, 

//i=fli(0;Q1,n2,Q3,fi4,a5)^fiA(0)1  (88) 

k= 1 

assuming  five  belts  of  sources.  Our  model  in  vector  notation  is  expressed  as 


dHmoi 

da 
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h  =  F(0)  +  s,  (89) 

where  h  ej$ixl  is  the  vector  of  scattered  magnetic  field  measurements  at  the  center  of  the  sensor  (L  is 
the  total  number  of  measurements  on  the  grid),  H  =  H(0)  is  the  Lx 5  matrix  as  a  function  of  0  and 

is  the  vector  of  normalized  magnetic  charges. 

Below  we  describe  three  gradient -based  methods  for  the  minimization  of  the  sum  of  squares. 

3.2.1  Stepwise  optimization 

This  is  a  traditional  way  to  recover  0  and  Q  in  alternating  fashion.  We  start  by  specifying  a 
stalling  value,  0() ,  and  compute  an  estimate  for  Q  from  linear  least  squares  by  solving  the  corresponding 
normal  equations, 

£2  =  (H'H)  ‘h'Ii,  (90) 

where  the  primes  denote  transposition.  We  then  keep  Q  fixed  at  Cl  and  perform  a  nonlinear  regression  to 
obtain  a  new  estimate  for  0  ,  which  we  then  use  to  find  a  new  Cl  using  (90),  iterating  until  we  attain 
convergence. 

3.2.2  Simultaneous  optimization 

In  simultaneous  optimization  we  treat  the  linear  Q  parameters  on  the  same  footing  as  the 
nonlinear  0  ;  the  Jacobian  thus  takes  the  form  of  an  Lx  10  matrix: 


5F 

5F 

5F 

5F 

5F 

5F 

5F 

5F 

5F 

5F 

5F 

5F 
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8Q 

dx 

5y 

dz 

d(p 

89 

an, 

5Q2 

5Q3 

5Q4 

5a 

3.2.3  Condensed  algorithm 

In  this  method,  we  use  a  closed-form  formula  to  reduce  the  number  of  model  parameters  from  10 

to  one: 

h  =  F(0)  +  s,  (92) 

where 

F(0)  =  h(h  h)  H'h  (93) 
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We  need  the  derivatives  of  F  with  respect  to  all  five  components  of  0  ;  to  compute  them  we  use  a  matrix 
chain  rule.  Using  a  formula  for  matrix  derivatives  we  get  the  derivative  for  a  specific  component  of  0  as 

—  =  — (H'H)  1  H'h  -H(H'H)  1  f — H  +  H—  ](H'H)  1  H'h  +  H(H'H)  1  — h  (94) 

<90  <90 V  ’  V  7  tv  <70  <70  y  ’  v  ’  <70 


Combining  the  five  derivatives  we  get  the  Lx 5  total  Jacobian 


<3F  GF  GF  SF  GF 
dx’  dy’  dz’  d(p'  36 


(95) 


and  the  Gauss-Newton  algorithm  takes  the  form 

es+1=es+(j;j  srlrsh  (96) 

to  be  iterated  until  reaching  convergence. 


3.3  Differential  evolution 

Differential  evolution  (DE)  [33-34],  one  of  the  global-search  algorithms  recently  developed  to 
bypass  the  local-minima  problem  that  often  leads  standard  gradient-search  approaches  to  make  incorrect 
predictions  for  location  and  orientation,  is  a  heuristic,  parallel,  direct-search  method  for  minimizing 
nonlinear  functions  of  continuous  variables.  Similar  in  concept  to  the  genetic  algorithms  that  have  been 
used  with  much  success  on  problems  with  discrete  variables,  DE  is  easy  to  implement  and  has  good 
convergence  properties. 

We  have  combined  the  DE  algorithm  with  the  above-discussed  dipole,  NSMS,  and  ONVMS 
techniques  to  invert  digital  geophysical  EMI  data  following  a  procedure  reminiscent  of  the  stepwise 
optimization  described  in  the  previous  section.  The  scattered  field  from  any  object  whose  location  and 
orientation  are  known  depends  linearly  on  the  magnitudes  of  its  responding  sources,  and  the  procedure 
starts  by  giving  initial  values  of  the  attitude  parameters  and  using  these  estimates,  along  with  the 
measured  data,  to  determine  the  source  amplitudes  by  solving  a  linear  system  of  equations.  The 
amplitudes  thus  found  are  fed  into  a  nonlinear  objective  function  that  quantifies  the  mismatch  between 
measured  data  and  model  predictions  and  whose  (DE-determined)  minimum  serves  to  refine  the  estimates 
for  location  and  orientation.  The  procedure  continues  to  alternate  between  these  linear  and  nonlinear 
stages  until  it  reaches  convergence  (or  a  preset  maximum  number  of  iterations).  The  responding 
amplitudes  are  then  stored  and  used  in  a  later  classification  step,  while  the  location  and  orientation 
parameters  are  used  during  target  excavation. 

Differential  evolution  uses  N p  -dimensional  parameter  vectors  v, 
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v pja'  P  1’^’  •••>  N p 


(97) 


where  G  is  a  generation/iteration  index.  In  our  case  v  =  {x  ,  yo,zQ,O,0} ;  the  first  three  are  the  object’s 
location  and  the  other  two  are  the  polar  ( 9 )  and  azimuthal  ( (f> )  Euler  angles  that  define  its  orientation  (by 
using  only  two  angles  we  are  assuming  that  UXO  are  effectively  BOR).  The  objective  function  to  be 
minimized  is  defined  as 


F(\)  = 


( MNfy 


M  N  f 

-Z2>, 

m=l  f=l 


f(v)-//dat;(v)  i2 

m,j  v  7  m,j  v  7 


(98) 


where  H*  (y)  and  /U:ij'  are  respectively  the  theoretical  prediction  (for  vector  v)  and  the  measured 
magnetic  field  data  at  the  m  -th  measurement  point  (of  M )  and  the  /  -th  frequency  or  time  point  (of 
N  f ).  The  DE  optimization  process  itself  can  be  subdivided  into  three  steps: 


1)  The  first  step  creates  random  initial  populations  v  ,  p  —  1,2,  ...,  N  ,  that  span  the  entire 
parameter  space.  For  a  given  \pC  in  the  generation,  a  linear  system  of  equations  is  constructed  by 
matching  measured  data  to  the  secondary  magnetic  field  from  (88).  This  system  is  linear  in  Q  and  is 
solved  directly  for  those  parameters. 

2)  The  second  step,  which  requires  the  most  execution  time,  is  the  calculation  of  the  secondary 
magnetic  field  (88)  for  each  of  the  v  pG .  When  the  NSMS  (or  ONVMS)  model  is  used,  the  calculation  for 

each  \pG  requires  a  fraction  of  the  time  required  to  execute  any  other  proposed  3D  forward  model;  this 

relative  computational  efficiency  makes  NSMS  (or  ONVMS)  an  attractive  alternative  for  performing  real¬ 
time  inversion. 

3)  Next  comes  the  evaluation  of  the  cost  function  for  each  population  member  and  the  storage  of 
the  best  sets  of  parameters.  At  each  step,  the  DE  algorithm  produces  an  estimate  of  position  and 
orientation.  By  examining  and  sorting  the  cost  function  at  each  step,  the  best-half  of  the  population  is 
chosen  as  the  next  generation’s  parameters,  whereas  the  bottom  half  is  discarded.  Thereafter  the  next 
generation  is  created  by  taking  the  parameters  in  the  previous  generation  and  applying  crossover  and 
mutation  operations  on  them.  The  three  steps  are  repeated  until  the  maximum  number  of  generations  has 
been  reached  or  until  the  objective  function  reaches  a  desired  value.  Rules  for  using  DE  are  discussed  in 
more  detail  elsewhere  [33-34], 
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3.4  The  HAP  method 


3.4.1  Estimating  the  location  and  orientation  of  buried  objects 

In  the  EMI  regime,  the  secondary  magnetic  fields  measured  by  the  EMI  receivers  are  induced  by 
eddy  currents  or  magnetic  dipoles  which  are  distributed  non-uniformly  inside  the  scatterer.  There  are 
some  particular  points,  named  “scattered  field  singularities”  (SFS),  where  most  of  these  sources  are 
concentrated.  Recent  studies  show  that  under  certain  conditions  the  entire  scatterer  can  be  replaced  with 
several  responding  elementary  sources  by  putting  them  at  SFS  points  [43-50].  The  mathematical  and 
physical  properties  of  SFS  and  its  applications  to  EM  scattering  problems  are  very  well  documented,  and 
their  study  is  known  in  the  literature  as  “ Catastrophe  Theory ”  [43-44].  Our  objective  has  been  to 
determine  the  locations  of  the  SFS  from  data  without  solving  traditional  ill-posed  inverse-scattering 
problems.  We  have  found  a  new  analytic  expression  for  estimating  the  location,  orientation,  and 
polarizability  elements  of  a  buried  object  stalling  from  measured  EMI  data.  The  algorithm  (dubbed 
“HAP”  [51])  is  based  on  the  fact  that  a  target's  response  can  be  approximated  by  dipole  sources 
concentrated  at  SFS  points.  It  utilizes  three  global  values  at  a  single  location  in  space:  (1)  the  magnetic 
field  vector  H,  (2)  the  vector  potential  A,  and  (3)  the  scalar  magnetic  potential  y/ .  Since  among  these 
quantities  only  the  H  field  (and  sometimes  only  one  of  its  components)  is  measurable,  we  employ  a 
variation  of  the  NSMS  model  to  obtain  A  and  t//  we  distribute  elementary  sources  on  an  auxiliary  planar 
layer,  located  between  the  sensor  and  the  object,  and  find  their  amplitudes  by  fitting  measured  data. 

The  magnetic  field  H  and  the  scalar  (  y/ )  and  vector  (A)  potentials  of  a  magnetic  dipole  are 


jkR 


H  = 


4  nR3 


f  3R(Rm)  ^ 


R 


m  j(l-jM)-k2(Rx(R  xm)) 


(99) 
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a  =  _a0 


m  x  R 

4  xR3 


(1  -jkRy* 


m  x  R 


G(  R) ,  where 


e'kR 

G(R)  =  — — (1-  jkR) 
4  n 


(100) 

(101) 


where  k  is  the  wave  number  in  the  surrounding  medium,  R  =  r  -  r  ,  r  is  an  observation  point,  and  r;  is 

the  location  of  the  dipole  [49]  (Figure  5).  Note  that  the  magnetic  field  (99)  has  terms  that  decay  as  R  l , 
R2 ,  and  R  ' .  The  range  kR  » 1  is  referred  to  as  the  far  zone,  and  fields  in  this  range  are  referred  to  as 
being  in  the  far  field.  Similarly,  fields  in  the  near  zone  kR  ^  1  are  referred  to  as  being  in  the  near  field, 
and  the  zone  kR&  lis  called  intermediate  zone.  Typically,  UXO  detection  and  discrimination  are 
conducted  in  the  near  zone.  In  addition,  in  the  EMI  regime  displacement  currents  are  considered 
irrelevant,  which  means  that  the  contribution  of  the  k2  term  in  equation  (99)  can  be  set  to  zero.  Making 
this  assumption,  taking  the  dot  product  of  (99)  with  R  ,  and  using  (100)  we  get  that 


H  R  =  H  (r-r,) 


R2  { 


R2 


(102) 


Similarly,  taking  the  cross  product  of  (99)  and  R  and  using  (101)  we  obtain 


HxR  =  G(«)T[3R(R,nll-nl)xR  =  -G(g,"lxR  A 


R'\  R 

Now,  the  cross  product  of  H  and  (103)  gives 


R  a0 


(103) 


Hx  — 
A0 


=Hx  [Hx  r]  =H(H  R)-R  |H|2  =2H  yr-R  |h|  , 


(104) 


which  allows  us  to  solve  for  R  : 

_2H  y/ - [H xA///# 


The  location  R  of  the  responding  dipole  is  seen  to  be  independent  of  the  frequency.  In  other 
words,  as  long  as  MQS  assumptions  hold,  equation  (105)  is  valid  when  the  dipole  is  in  free  space  and 
equally  well  when  it  is  embedded  in  a  conducting  medium  such  as  seawater.  Also  note  that  R  is 
determined  as  a  ratio,  which  makes  the  expression  (105)  partially  tolerant  to  noise  due  to  scaling 
arguments,  since  A  and  t//  are  dependent  on  the  H  field  (see  equations  (102)  and  (103)).  Taking  the 


(105) 
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cross  product  of  R  and  (103)  from  the  left  side  and  using  equation  (102)  we  obtain  an  expression  for  the 
dipole  moment  m : 

m"G^)(R^  +  fA///,,XR5  (1°6) 

with  R  previously  determined  from  equation  (105). 

3.4.2  A  simplified  HAP  method 

It  is  possible  to  simplify  the  HAP  method  by  eliminating  the  need  for  the  vector  potential.  We 
rewrite  equation  (102)  as 

Hrrf  =-2^+H-r,  (107) 

which  provides  a  least-squares  estimate  of  rd  when  evaluated  at  N  distinct  observation  points: 

H  (r.)  H  (r.)  H  (r.) 

xx  1 7  y  v  1  7  z  v  1 7 

H  (r,)  H  (r  )  H  (r7) 

2;  y  v  2  7  z  v  2  7 

H  (r„)  H  (r„)  H  (rM) 

xy  N  7  yx  N  7  zx  N  7 

3.4.3  Determining  the  HAP  amplitudes 

To  construct  the  potentials  (and  the  other  field  components,  if  unavailable)  we  assume  that  the 
field  is  produced  by  a  surface  distribution  of  magnetic  charge  q(s')  spread  on  a  fictitious  plane  located 

just  below  the  ground  (Figure  6).  The  positions  r  of  the  sources  are  fixed  and  known  by  construction, 
and  the  field  can  be  expressed  as  the  matrix-vector  product 

Hz(r)  =  \q{s'\Z-Z’'  ds'  =  Zz-q  (109) 

J  i*  —  r 

by  employing  a  quadrature  scheme.  To  determine  the  array  q  of  charges  we  minimize  the  difference 
between  model  predictions  and  collected  data  Hmeas  at  a  set  of  known  points: 

q  =  argminj(Z;.q-H— )2=[Z^Z:]_1[Z^H— ],  (110) 

where  each  matrix  row  corresponds  to  a  different  measurement  point  and  each  column  to  a  subsurface  of 
the  underground  virtual  source  layer.  The  potential  is  then  found  from 
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<nl) 

Current  EMI  sensors  operate  in  both  monostatic  and  multistatic  modes.  Monostatic  sensors,  such 
as  the  Geophex  frequency-domain  GEM-3  instrument  [48]  and  the  Geonics  EM-61  and  EM-63  time- 
domain  instruments  [47]  have  collocated  transmitter  and  receiver  coils,  whereas  multistatic  sensors  like 
the  MPV  time-domain  instrument  [28]  and  the  Berkeley  UXO  Discriminator  (BUD)  [52]  have  multiple 
transmitters  or  multiple  receiver  coils  or  both.  We  have  implemented  numerical  procedures  to  estimate 
the  vector  and  scalar  magnetic  potentials  starting  from  multi-static  or  mono-static  EMI  data.  For  bistatic 
data  we  determine  the  potentials  as  described  above;  for  the  monostatic  case  we  normalize  the  amplitudes 
of  the  responding  auxiliary  sources  by  the  primary  magnetic  field.  The  procedure  is  discussed  in  further 
detail  in  [51]. 

It  is  worth  reiterating  that  the  HAP  method  replaces  the  scatterer  with  a  point  dipole,  and  is  thus 
based  on  a  rather  drastic  simplification;  yet  it  provides  acceptable  location  estimates  because  the  sources 
within  the  target  that  produce  the  scattered  field  tend  to  concentrate  at  a  set  of  “scattered  field 
singularities”  [46,  50] .  The  locations  of  these  singularities  change  at  every  measurement  point,  since  the 
primary  field  of  the  sensor  also  changes;  the  HAP  method  takes  these  variations  into  account  and  outputs 
an  average  location  as  a  result. 


Figure  6:  Determining  the  location  and  orientation  of  a  buried  target.  The  method  assumes  the  object  is  a  point 
dipole  and  exploits  an  analytic  relation  between  the  field  measured  at  r  and  the  scalar  potential  at  the  same  point  to 

find  the  location  rd .  The  potential  is  constructed  using  a  layer  of  equivalent  magnetic  sources  placed  between  the 
sensor  and  the  object;  r  is  a  typical  location  on  the  layer. 
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3.4.4  The  HAP  method  with  gradient  information 

The  HAP  technique  can  be  simplified  further  by  reducing  the  formulation  such  that  it  only 
requires  the  magnetic  field  and  its  gradient,  both  of  which  are  measurable  by  current  sensors.  After  taking 
the  gradient  of  equation  (102)  with  respect  to  the  x-,  y-,  and  ^-coordinates,  we  obtain 


BHx  8H ,  dH .  8Hx  oH  dH 

xd^-  +  yd^-  +  zd—^  =  3Hi+x  —  +  } ’  —  +z—^- 
ox  ox  ox  ox  ox  ox 

8HK  8H  dH.  GH  8H  dH 

8y  8y  8y  8y  8y  8y 

8H  8H  dH.  8H  8Hy  dH. 

xd^  +  yd^-  +  zd—^  =  3Hz+x—^  +  y  —  +z—^- 
l  6z  8z  8z  8z  8z  8z 


Thus,  in  order  to  determine  the  target’s  location  we  need  only  the  magnetic  field  H  and  its 
gradient  at  a  given  point  in  space. 
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4  Statistical  approaches 

4.1  Introduction 

UXO  discrimination  is  a  three-step  process  that  comprises  data  collection,  parameter  inversion, 
and  target  classification.  The  classification  step  requires  robust  statistical  approaches  that  will  dependably 
separate  UXO  from  innocuous  items.  As  inputs  to  these  procedures  one  needs  inverted  intrinsic 
parameters  that  encapsulate  relevant  features  of  the  targets  in  question  (like  the  time  decay  of  their 
response,  any  symmetries  they  may  possess,  etc.)  in  such  a  way  that  the  parameters  remain  as  consistent 
as  possible  between  samples  of  the  same  type  of  target  and  vary  significantly  enough  from  type  to  type.  In 
other  words,  one  needs  to  find  parameters  that  cluster  well.  One  can  then  create  feature  libraries  and, 
given  an  unknown  target,  compare  its  features  to  those  in  the  library  and  ascribe  it  to  the  class  or  cluster 
that  better  assimilates  it.  There  are  several  clustering  techniques  available,  such  as  K-means  clustering 
[53],  principal-component  analysis  (PCA)  [21-25,  54],  or  the  support  vector  machine  (SVM)  algorithm 
[55-56],  that  are  largely  heuristically  motivated  and  do  not  require  an  underlying  statistical  model.  A 
possible  alternative  is  model-based  clustering,  which  assumes  that  the  samples  of  each  type  of  target  in  an 
inverted  feature  data  set  follow  identical  multivariate  statistical  distributions  (Gaussians,  for  example)  but 
with  different  parameters,  and  that  the  complete  set  is  a  mixture  of  a  finite  number  of  such  distributions. 
This  model -based  approach  has  the  desirable  traits  (l)that  it  permits  the  use  of  objective  statistical 
criteria — like  the  Akaike  Information  Criterion  (AIC)  or  the  Bayesian  Information  Criterion  (BIC) — to 
choose  the  optimal  number  of  clusters  and  (2)  that  there  are  several  available  distributions  that  can  be 
varied  until  one  is  found  that  best  fits  the  data;  for  “heuristic”  algorithms,  on  the  other  hand,  choosing  the 
“correct”  number  of  clusters  and  the  best  clustering  method  is  not  a  settled  question  and  has  to  be  dealt 
with  on  an  ad  hoc  basis. 

Clustering  methods  can  be  either  unsupervised  or  supervised.  Supervised  clustering  uses 
parameters  derived  from  training/calibration  samples,  whereas  unsupervised  clustering  applies  the  same 
classification  criteria  to  all  targets,  regardless  of  size,  composition,  and  decay  curves.  Supervised 
clustering  has  the  obvious  advantage  that  it  utilizes  additional  information  from  the  training  data  as  prior 
knowledge  and  the  obvious  drawback  that  it  completely  ignores  potentially  valuable  information  from  the 
blind  data.  Several  supervised  clustering  techniques  (like  SVM  [57-58]  and  template  matching)  have  been 
applied  to  UXO  discrimination.  A  straightforward  approach  to  implement  the  model-based  supervised 
clustering  algorithm  is  to  (a)  estimate  the  parameters  from  the  Paining  samples  and  (b)  use  the  estimated 
values  of  the  parameters  to  classify  blind  data. 
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This  chapter  introduces  and  analyzes  several  such  methods  from  a  UXO  discrimination  point  of 
view,  including  mixed  modeling,  supervised  and  unsupervised  Gaussian  mixture  models,  and  the  SVM 
technique.  We  first  introduce  and  discuss  mixed  modeling,  then  give  an  overview  of  the  unsupervised 
multivariate  normal  mixture  model,  and  end  by  presenting  the  SVM  approach  to  target  classification. 


4.2  Mixed  modeling  applied  to  advanced  EMI  features 

Mixed  modeling  may  be  viewed  as  a  compromise  between  a  frequentist  and  a  Bayesian  approach. 
As  in  a  Bayesian  approach,  it  formulates  the  problem  in  a  hierarchical  fashion.  The  first  stage  is  a 
conditional  model  that  describes  the  measurements  as  a  function  of  the  signal,  assuming  that  the 
distribution  of  the  true  signal  is  known.  The  second  stage  specifies  the  distribution  of  the  true  signal  using 
a  distribution  given  a  priori.  In  mixed  modeling,  as  in  the  frequentist  approach,  the  a  priori  distribution  is 
given  but  the  parameters  are  unknown.  These  parameters  are  further  estimated  using  maximum  likelihood 
or  an  approximation,  especially  when  the  dimension  is  high,  to  avoid  ill-posedness.  Under  this  project  we 
investigated  (1)  the  applicability  of  mixed  modeling  to  regularize  the  inversion  of  EMI  data  and  (2)  the 
possibility  of  expanding  the  procedure  to  incorporate  target  discrimination. 

In  general,  the  determination  of  intrinsic  features  of  targets  (dipole  polarizabilities,  NSMS, 
ONVMS,  or  the  like)  is  an  ill-posed  problem  that  often  results  in  instability  of  inversion.  One  may 
improve  this  instability  by  incorporating  a  priori  information  about  the  properties  of  the  relevant  objects. 
The  method  we  suggest  here  uses  the  total  NSMS  of  Section  2.3  for  illustration  purposes;  the  total 
ONVMS  of  Section  2.4  works  just  as  well.  We  start  by  writing  down  an  expression  for  the  secondary 
field  due  to  a  set  of  discrete  sources 


H 


'  (ry  )=Z n,G(r, . : rj)  .  where 


G(r „r.)  =  Hpr(r  ,r.)  —  J  -?-(3r..(r..  -n  )-r2n.)  j 

V  t  J'  n  v  o’  i'  ^  5  V  Ji V  Ji  i  7  Ji  ' 7  1 


(113) 


is  the  Green  function,  rfl,  r y,  and  r,  are  respectively  the  locations  of  the  o-th  transmitter,  the  j-th  receiver, 
and  the  i-th  NSMS  source,  r  =  r  -  r  ,  and  Ns  is  the  total  number  of  NSMS  sources.  The  total  NSMS  (or 
ONVMS)  is  the  main  classification  parameters;  we  will  focus  directly  on  estimating  the  total  NSMS  using 
the  average  amplitude  Q .  Mixed  modeling  relates  the  average  Q  to  the  actual  amplitude  at  each  k- th 
point  through 

Q,=0  +  £,.,  (114) 


where  the  si  are  random  deviations  from  the  average  with  mean  £(£\)  =  0  and  variance  varU;  )  =  a]: . 
After  combining  equations  (113)  and  (114),  and  matching  modeled  field  to  the  data  we  have  that 
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dj  ■  £)^c;(r.r, )  +  //.. 


(115) 


where  dj  is  the  data  collected  at  the  j-th  point,  j  =  1,  2,  3,  Nd,  and  //  is  a  random  variable  with  zero 
mean  and  variance 

var(//; )  =  cx 2  +  cr  G2  (r. ,  r . ) .  (116) 


where  cr2  is  a  parameter  that  is  determined  using  the  mixed  model  approach.  This  is  the 
fundamental  principle  of  mixed  modeling,  which  assumes  that,  in  the  same  way  that  data  can  be 
contaminated  with  random  noise,  the  model  used  to  process  and  invert  it  can  also  contain  random 
statistical  errors.  The  error  statistics  for  both  data  and  model  are  estimated  using  the  mixed  model,  which 
belongs  to  the  family  of  nonlinear  marginal  models  (Chapter  6,  [59]).  The  average  Q.  can  be  estimated 
from  actual  data  dj  using  a  regularized  nonlinear  least  squares  minimization: 

f  ns  ^ 2  n, 

^\drt^G(  r.,r.)  +  p£(Q.-Q)2.  (117) 

i,j  V  i  J  (=1 


The  regularization  parameter  p  is  found  by  regressing  the  squared  residuals  on  the  variance  var(  //  )  as  a 
linear  function  of  <jg  .  The  magnitudes  Q  of  the  NSMS  are  determined  by  solving 


Q  =  (pi  +  H'H)_1(pQl  +  H'y) . 


(118) 


Note  that  we  use  matrix  notation  (boldface  type  for  vectors  and  matrices,  regular  roman  type  for  scalars, 
and  primes  to  denote  transposition).  The  mixed  model  (MM)  algorithm  for  finding  Q  can  be  summarized 
as  follows: 


1 .  Estimate  the  NSMS  from  (118)  using  standard  nonlinear  least  squares  and  taking  p  —  0  . 

1  N'; 

2.  Obtain  Q  from  Q  =  —  V  Q  . 

A hi 


Ns 

3.  Obtain  the  squared  residuals  e2  and  regress  them  on  y  Gj(r,r  )  to  determine  estimates  of 


cj;  and  cr: .  Let  p -a]  lap 

4.  Estimate  the  NSMS  from  (118)  and  repeat  steps  2,  3,  and  4  until  Q  converges. 
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We  now  proceed  to  describe  an  algorithm  to  obtain  estimates  of  Q ,  cr2 ,  and  al  using  a 
maxi  mum-likelihood  approach.  We  rewrite  the  NSMS  model  (115)  as: 

d  =  GQ  +  ?7,  Q  =  Ql  +  s,  (119) 

where  d  =  (dl,d2,...,dN  )'  is  the  vector  of  all  Nd  measurements,  G  gRNjXN°  is  a  matrix,  il  and 

1  =  (1,1,...,  N  )'  are  vectors  of  size  N  x  1 ,  and  the  average  amplitude  Q  is  a  scalar.  We  assume  that  the 
errors  are  normally  distributed: 

y  -  m n  h, ,  a2£lNd  +  ^2GG') ,  (120) 

where  N  denotes  a  multivariate  normal  distribution  and 

hj  =  G1 ,  (121) 

which  implies  that 

E(y)  =  Q  hr  and  cov(d)  =  erl  +  crGG' .  (122) 

In  previous  studies  we  used  weighted  least  squares  and  estimated  the  regularization  parameter 
using  the  method  of  moments.  Here  we  apply  the  method  of  maximum  likelihood,  which  is  more  complex 
but  more  precise.  To  obtain  the  log-likelihood  we  use  the  following  fact  of  the  multivariate  statistics:  if 
d  ~  /V(p,  fi) ,  then  the  density  function  is  [60] 

p(d;p,a)  =  (2^r/2|nf  1/2  exp  -|(d-n)'n-1(d-(i)  .  (123) 

The  likelihood  function  is  the  density  function  with  data  and  parameters  exchanged, 

L(n,«)  =  ( 2^)“,,/2|Qf1/2exp  -|(d-n)'n_1(d-rt  (124) 


The  twice-negative  log -likelihood  (TNLL)  function  -2  In  L  takes  the  form 

In  |cr2I^  +  <t2GG'|  +  (d  -  Qh;  )'(cr2 1^  +  cx2GG')^  (d  -  Qh, ) .  (125) 

2  2 

Minimizing  the  TNLL  yields  maximum-likelihood  estimates  of  a e  and  crs  ,  and  once  these  are  known 
we  obtain 


n  = 


d'(*X  +*s  GG')-'h, 
h[(cr2X  +<t2sGG'TX 


(126) 
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The  straightforward  minimization  of  this  problem  may  be  prohibitive  because  matrix  Iw  +  vGG'  has 
dimension  N  x  N  d,  so  we  use  a  reduced-dimension  matrix  inverse  [59] 

(I*,  +  vGGT1  =  1^  -  G(/oI  v  +  G'G)  G' ,  (133) 

where  p- cr  /  <j2g  = 1/ v  is  the  regularization  parameter  [61-65].  Thus  we  do  not  have  to  invert  a 
N  d  x  N  d  matrix  but  only  one  of  size  NsxNs.  Similarly,  we  can  reduce  the  computation  of  Eq.  (132): 

In  I  +  vGG'  =  In  plN  +  G'G  -  Ns  In  p .  (134) 

Finally  we  come  to  the  one-dimensional  minimization  (one -variable  function): 
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f  ip)  =  Nd  ln(||r||2  -g'(pIWj  +  G'G)-1  g)  +  In  |pl  ^  +  G'G| -7VS  lnp 
where  g  =  H'r  for  p>  0  .  Equivalently  one  can  solve  the  equation 


(135) 


df 


N, 


N 


J  „  l|2  - £'(pIN  +G'G)“2g  +  tr(/7lw  +G'G)_1 - *  =  0, 

dp  r'-gVIv  +G'G)-1g  '  1  P 


(136) 


which  can  be  expressed  as  H(p)  =  0,  where 


H(p)  =  Ndpg'(pl„  +  G'G)“2g+  ptv(p\N  +  G'G)"1  -  TV, 


r  ~-g'(jA  +G'G)_lg 


(137) 


But 


^  -  ng'(plN  +  G'G)-2  g  -  2npg'(pl  N  +  GG)3g 

dp 


[tr (plNi  +G'G)1  -ptr(plNs  +GG)-2][|rf  -g'(plNt  +G'G)  'g 
[ptr(plNl  +G'G)-1-^]g'(pINj  +GG)-2g 


(138) 


Therefore,  the  Newton’s  gradient  search  algorithm  takes  the  form 

H(PS) 


Ps+i  Ps 


dH  /  dp 


(139) 


starting  from  pQ  =  0  .  Finally,  the  algorithm  is  as  follows: 

1 .  Obtain  the  least-squares  Q  as  y'h  /  hjh  . 

2.  Compute  r  =  d  —  Qht  and  g  =  G'r . 

3.  Find  the  minimum  of  f(p)  using  Newton’s  p -iterations  above. 

4.  Compute  the  new  Q  as 


d' 

I  -G(/Jl„ +G'G)  G^ 

h,  d'h-d'G(/2lw +G'G)  ^G'h, 

h; 

ln-G{p\N +G'G)G' 

h[  ”h;h,-h;G(pIw  +G'G)1G'h| 

and  go  to  step  2.  Incidentally,  this  is  the  same  formula  (126)  but  expressed  in  terms  of 
p  -  cr2  /  cr2 .  Indeed,  we  have: 
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_  _  d'(I^  +  <f2g  /  ^GG')-1h1  _  d'(I^  +  vGG  )  h, 
Q”h;(I^+^/^GG')-|h|  ■h;(IiVj  +  vGGTIh1  ■ 

5.  Iterate  until  convergence,  and  at  the  end  compute  the  amplitudes: 

£2  =  {plN  +  G'G)  1  (G'd  +  pQl) . 

4.3  Gaussian  mixture  models 

4.3.1  Model-based  supervised  clustering 

Targets  of  interest  (TOI)  with  similar  features  (i.e.,  dipole  polarizabilities  or  total  NSMS  or 
ONVMS)  are  likely  to  show  similar  power-law/exponential  time  decay  patterns  under  various  conditions, 
and  as  a  result  these  patterns  form  clusters  when  plotted  in  a  convenient  and  pertinent  space.  It  is  possible 
to  identify  an  unknown  target  by  comparing  its  time -dec  ay  parameters  to  those  of  a  set  of  previously 
characterized,  previously  clustered  objects  and  assigning  it  to  the  category  where  its  profile  fits  best.  Such 
“supervised”  clustering  allows  the  use  of  additional  information  from  the  training  data  as  prior 
knowledge;  on  the  other  hand,  it  uses  only  the  training  set  to  estimate  the  parameters  and  completely 
ignores  the  blind  data,  which  is  potentially  quite  useful.  The  test  data  set,  moreover,  is  usually  much 
larger  than  the  training  sample,  implying  that  the  unused  information  may  be  substantially  richer  than 
what  is  contained  in  the  training  sample. 

Let  us  assume  that  there  are  K  clusters  and  that  each  cluster  is  mathematically  described  by  a 
parametric  continuous  or  discrete  distribution  function  (usually  a  Gaussian,  as  in  this  case).  The 
classification  parameters  (extracted  for  example  by  fitting  the  total  NSMS  or  ONVMS  with  a  Pasion- 
Oldenburg  time-decay  law)  can  then  be  arranged  in  an  n  x  m  matrix  denoted  by  Y  =  [Y|,  Y2,  ...,  Ym], 
where  Y„  i  —  1,  2,  ...,  n,  is  a  vector,  n  is  the  number  of  anomalies,  and  in  is  the  number  of  parameters. 
Each  Y,  can  be  considered  to  follow  an  m-dimensional  mixture  of  normal  distributions  expressed  as 

(143) 

k= 1 

where  wk  is  the  mixing  weight  of  cluster  k  (defined  as  the  proportion  of  anomalies  that  belong  to  it), 

Ki wt  = 1  >  and 

f,(Y,  1 IV  )  =  I  1  exp[““ (Y.  - PA, I'd’1! Y.  - )  I  (144) 

W2  *)"  1  2  2 


(141) 


(142) 


Sky  Research,  Inc. 


48 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


is  the  probability  density  of  the  k- th  normal  distribution  with  the  in  x  1  mean  vector  p.  and  the  m  x  m 
variance-covariance  matrix  a.  .  The  parameters  p  ,  ct;  ,  and  wk  are  estimated  by  the  maximum 
likelihood  (ML)  criterion  using  the  Expectation  Maximization  (EM)  algorithm. 

This  simple  and  intuitive  supervised  method  usually  performs  well  if  the  number  of  TOI  within 
each  known  cluster  in  the  training  sample  is  sufficiently  large  to  ensure  high  accuracy  of  the  estimates  of 
pt  and  a .  .  For  small  Paining  samples  these  estimates  are  subject  to  large  errors — in  some  cases,  if  the 

number  of  targets  within  a  cluster  is  smaller  than  the  number  of  parameters,  the  estimated  variance- 
covariance  matrices  may  not  even  be  positive  definite.  Furthermore,  much  information  from  the  test 
dataset  has  not  been  fully  utilized.  The  test  dataset  is  usually  much  larger  than  the  Paining  sample, 
implying  that  the  unutilized  information  may  be  substantially  more  than  that  contained  in  the  training 
sample.  To  overcome  this  problem  we  use  unsupervised  clustering  and  derive  p  and  cr  from  blind-test 
data  using  an  iterative  EM  algorithm. 

4.3.2  Unsupervised  classification  using  the  multivariate  normal  mixture  approach 

Mixture  distribution  is  perhaps  the  only  model-based  approach  among  existing  methods  of 
clusterization  and  pattern  recognition.  Its  attractive  features  are  that  (a)  it  is  not  necessary  to  specify  what 
class  each  observation  belongs  to  (i.e.,  the  classification  is  “unsupervised”)  and  that  (b)  the  method 
estimates  the  membership  probability  which  results  in  confusion  matrix  n  . 

Let  x  ,...,x  be  m-dimensional  feature  vectors  that  we  want  to  split  into  K  classes.  It  is  assumed 
that  each  x  belongs  to  one  of  K  classes  that  are  described  by  densities  <p  (x),...,<p  (x) .  If  n  >0 
denotes  the  probability  of  x  belonging  to  class  k  ,  then  the  mixture  density  is  the  linear  combination 

K 

(p{x)  =  YJttk(pk{x),  (145) 

k= 1 

where 

I>*=  L  (146) 

k=l 

In  the  case  of  a  normal  distribution  we  have 

f>,ATx;p,,Q,)  (147) 

k= 1 
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for  the  mixture,  where  the  mean  and  the  covariance  matrix  £lk  are  different  for  the  different  clusters 
and  subject  to  estimation  along  with  probabilities  n .  Since  the  density  of  ATx;p,  ,Q,  )  is 


%(x;p,,n,)  =  (2^)-'"/2|n,|-1/2exp 


]_ 

2 


(x-nt)'n^(x-nt) 


(148) 


the  method  of  maximum  likelihood  prescribes  the  maximization  of  a  nonlinear  function: 

L(HkSlk,7Tk)  =  (2 nTmllY]tk  In*r1/2  exP 

k=\ 

This  task  is  not  easy  because  the  function  is  unbounded  (it  may  go  to  +co ),  and  this  creates  numerical 
obstacles. 

A  penalized  version  of  the  multivariate  normal  mixture  method  was  recently  developed  [66] .  The 
idea  of  penalization  is  based  on  the  fact  that  the  problem  of  mixture  maximization  is  not  difficult  when  all 
Clk  are  the  same,  and  thus  it  would  make  sense  to  penalize  for  variation  among  the  Q  .  Moreover,  when 

C\  -  Cl  then  they  should  be  equal  to  the  sample  covariance  matrix 

A  =  Sr  =  -£(x,.  -x)(x;.  -x)' .  (150) 

n  ;=1 

Specifically  Chen  and  Tan  suggest  to  penalize  the  log  of  (148)  with: 

P  =  -“,E[«S,‘/.;1)  +  ln|'/.,|].  (151) 

k=\ 

The  justification  of  this  penalization  is  that  it  reaches  an  extremum  if  and  only  if  C\  =  S  so  its 
addition  to  (148)  tends  to  keep  the  Clk  close  to  S  .  The  authors  suggest  using  either  an=l/n  or 

an-H  4n  for  the  weight  coefficients.  Once  penalization  is  in  place  we  apply  the  expectation- 
maximization  (EM)  algorithm  to  find  the  maximum  of 

l(\lk,nk,7t)  + p  (152) 


^(x-^y£V(x-^) 


max 

-Kk 


(149) 


over  Kk,  PrA’  where  the  log-likelihood  /  =  In  L .  Specifically,  the  EM  algorithm  for  the  multivariate 
normal  mixture  (MNM)  is  as  follows: 

1.  Let  the  initial  estimates  of  7rf , p , , Q  be  given  (for  example,  we  can  take  /r  =  1  IK, 
Clk  -  S  t ,  and  p;  from  K ). 
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2.  Compute  the  elements  of  the  confusion  matrix  n  as  the  probability  that  a  given  observation 
x  belongs  to  class  k  : 


_ 7C<p,(x;p,,f2,) 

%  -  - - 

^?.(x;n.,n.) 

j= i 

3.  Adjust  the  probabilities  accordingly: 


4.  Recompute  the  means 


and  covariance  matrices 


where: 


n. 


K 


ik 


TVx 

4-^  ik  i 


i=l _ 

TUT, 

k 


_  2fl„S,-+St 

2  an  +  nnk  ’ 


i=i 


(153) 


(154) 


(155) 


(156) 


(157) 


5.  If  the  values  are  different  from  the  previous  iteration  return  to  Step  1  and  continue  until 
reaching  convergence. 


4.4  Support  vector  machines  for  subsurface  object  classification 

The  Support  Vector  Machine  (SVM)  [38]  is  a  machine-learning  methodology  based  on  statistical 
learning  theory  that  has  been  used  to  perform  binary  classification  [67]  and  regression  [68]  and  has 
recently  been  adapted  for  multi-category  classification  [69].  The  method  has  been  employed  in  UXO 
research,  either  to  classify  or  regress,  in  combination  with  the  point -dipole  model  [7,  70-71],  the 
Standardized  Excitation  Approach  [56,  72-73],  and  finite  elements  [74-75],  and  has  shown  to  be 
competitive  in  its  discrimination  ability  in  relation  to  neural  networks  [73-74]  and  other  statistical 
methods  [57,  76]. 
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A  support  vector  machine  learns  from  data:  when  fed  a  series  of  answered  training  examples  it 
attempts  to  make  sense  of  them  by  weighing  the  available  empirical  evidence,  with  no  need  for  an 
underlying  model,  and  applies  this  knowledge  to  make  predictions  about  unseen  cases.  The  examples  can 
be  any  combination  of  model  parameters  expected  to  contain  evidence  of  the  essence  of  an  object.  In  the 
simplest  instance  of  binary  classification  each  n-dimensional  example  x  has  an  associated  yes/no 

attribute  y.  —  ±1 ;  the  SVM  performs  the  classification  by  finding  a  hyperplane  that  divides  the  parameter 

space  into  two  distinct  regions,  each  of  which  ideally  contains  points  from  only  one  of  the  categories. 
During  the  learning  or  training  process  the  machine  readjusts  the  hyperplane  parameters  to  accommodate 
every  training  vector  until  it  strikes  an  optimal  balance  between  fitting  accuracy  and  model  simplicity.  All 
information  about  the  hyperplane  is  contained  in  a  subset  of  the  examples — the  support  vectors  that  give 
the  method  its  name — which  are  then  combined  to  specify  a  predicting  function. 

The  SVM  algorithm  uses  two  different  strategies  to  tackle  the  nonseparability  of  realistic  data.  On 
the  one  hand  it  projects  the  examples  into  a  space  of  higher  dimensionality  by  means  of  a  kernel  function 
[77].  The  separating  surface  thus  found  is  flat  by  construction  in  the  new  space  but  can  be  curved  and 
even  multiply  connected  in  the  original.  On  the  other  hand,  the  technique  hies  to  control  overfitting — and 
thus  concentrate  on  essentials  rather  than  on  details,  resulting  in  better  generalization — by  having  an 
adjustable  penalty  on  misclassifications.  This  penalty  is  represented  by  a  single  scalar  parameter,  the 
capacity  of  the  machine  [78]. 

During  training  an  SVM  solves  the  constrained  quadratic  optimization  problem  [79] 


max 

a 

S.t. 


i  ^  i,j 

=  0,  0  <  a,  <  C, 


(158) 


whose  solution  is  a  vector  of  coefficients  a  that  measure  the  information  content  of  the  examples  and  are 

nonzero  only  for  the  support  vectors.  The  coefficients  are  prescribed  not  to  exceed  the  capacity  C  ,  which 
limits  the  influence  of  potentially  problematic  points  on  the  final  result. 

The  projection  to  higher  dimensions  occurs  by  substituting  the  scalar  products: 

xjxj  ->/s:(xi.,x.)  =  /(xi)T/(x.)  (159) 


for  some  mapping  <j)(x) .  The  function  K  is  the  kernel  we  mentioned  earlier.  It  is  not  necessary  to  know 
<|)(x)  to  find  K  :  any  function  that  combines  two  vectors  into  a  scalar  and  fulfills  the  (not  very  restrictive) 
set  of  conditions  spelled  out  in  Mercer’s  theorem  [80]  can  be  used  as  a  kernel.  Some  kernels  stretch  out 
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the  examples  into  the  added  dimensions  in  such  a  way  that  gaps  open  up  between  the  examples  which 
permit  a  flat  separating  surface  to  pass  through.  In  this  chapter  we  use  the  radial  basis  function  (RBF) 
kernel: 


K(x.,x.)  =  exp(— (x.  -x  ,)T(x.  -xp/2cr2), 


(160) 


which  surrounds  every  example  with  a  surface  that  in  a  sense  “repels”  the  separating  hyperplane.  The 
Gaussian  width  cr  is  a  second  adjustable  parameter  and  usually  has  a  scale  on  the  order  of  the  average 
separation  between  points.  In  [81]  it  was  found  that  polynomial  kernels  may  outperform  the  RBF  kernel 
in  some  electromagnetic  inverse  problems.  We  find  that  the  linear  kernel  makes  similar  predictions  and 
runs  faster  than  the  RBF,  though  the  difference  in  run  time  is  negligible  for  the  number  of  training  data 
and  example  features  that  we  typically  use.  Once  a  is  known  the  SVM  can  predict  the  attribute  of  an 
unknown  example  using  the  function  [78,  82] 


/ (x) = sgn 


Eaji^(X/,x) 

VieSV  J 


(161) 


There  are  several  ways  to  generalize  the  SVM  procedure  to  perform  multi -class  categorization. 
These  have  been  reviewed  in  [69],  whose  authors  conclude  that  the  methods  more  suitable  for  practical 
use  perform  several  binary  classifications  instead  of  attempting  to  separate  all  classes  at  once.  In  this 


work  we  adopt  a  one-against-one  approach  [82]  in  which  the  system  carries  out 


k(k- 1) 
2 


optimizations  and  obtains  the  same  number  of  decision  functions  of  the  form  (161).  When  given  an 
example  to  predict  the  algorithm  proceeds  by  ballot:  it  evaluates  the  decision  functions  one  by  one  on  the 
example  and  adds  a  vote  to  the  one  category  (out  of  two)  in  which  it  is  predicted  to  be.  At  the  end  the 
example  is  assigned  to  the  category  with  the  most  votes;  should  there  be  a  tie  between  two  classes,  the 
program  arbitrarily  selects  that  with  the  smallest  label. 
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5  Advanced  models  applied  to  next-generation  sensors: 

Modeling  and  validation 

5.1  Introduction 

A  wide  range  of  different  electromagnetic  induction  sensing  technologies,  with  novel  waveforms, 
multi-axis  transmitters,  and  scalar/vector  receivers  have  been  recently  developed  under  SERDP-ESTCP 
programs.  These  advanced  EMI  sensors — including  the  MetalMapper,  the  TEMTADS  array,  the  Berkeley 
UXO  discriminator  (BUD),  and  the  man-portable  vector  (MPV)  sensor — provide  measurements  that 
feature  a  combination  of  high  spatial  diversity,  different  viewpoints,  and  a  very  wide  dynamic  range  and 
which  do  full  justice  to  the  vector  character  of  the  electromagnetic  field.  Current  state-of-the-art  EMI 
systems  thus  offer  data  of  unprecedented  richness  for  use  by  discrimination  processing  algorithms.  We 
have  adapted  our  advanced  EMI  models  and  data-interpretation  and  -processing  schemes  to  all  these 
innovative  EMI  systems  in  order  to  take  advantage  of  the  quality  of  the  data  they  provide. 

This  chapter  overviews  these  advanced  EMI  sensors,  their  geometries  and  sensing  modalities,  and 
the  procedures  we  have  in  place  to  model  the  way  they  establish  primary  fields  and  measure  subsurface 
responses.  We  validate  our  methods  by  making  comparisons  between  measured  and  modeled  data  for 
single-  and  multi-target  scenarios.  We  initially  describe  the  MetalMapper,  continue  with  TEMTADS  and 
BUD,  and  finish  with  a  look  at  the  MPV. 

5.2  MetalMapper 

The  MetalMapper  (MM)  is  an  advanced  EMI  system  for  UXO  detection  and  discrimination 
developed  primarily  by  G&G  Sciences  and  commercialized  by  Geometries.  The  system  has  three 
mutually  orthogonal  transmitter  rectangular  loops.  It  is  able  to  illuminate  a  target  with  primary  fields  from 
three  independent  directions  from  a  single  spatial  field  point.  The  1  m  x  1  m  Z  transmitter  loop  is  located 
at  ground  level.  The  Y  transmitter  loop,  also  1  m  x  1  m,  is  centered  56  cm  above  the  Z  loop,  as  is  the 
0.98  m  x  0.98  m  X  transmitter  (Figure  7).  The  targets  are  illuminated  from  different  directions  depending 
on  the  geometry  between  a  particular  transmitting  loop  and  the  target.  The  system  has  seven  10-cm-side 
receiver  cubes  placed  at  seven  unique  spatial  points  on  the  plane  of  the  Z  transmitter  loop.  The  receivers 
measure  the  vector  cIB  I  dt  at  each  of  the  seven  points,  thus  providing  63  independent  readings  of  the 
transient  secondary  magnetic  field  for  each  instrument  location.  The  positions  of  the  receiver  cubes’ 
centers  with  respect  to  the  Z  transmitter  loop  (whose  center  we  consider  as  the  local  origin  of  coordinates 
for  the  system)  are  given  in  Table  1. 
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Figure  7:  The  MetalMapper  during  SLO  site  deployment  (left)  and  its  schematic  diagram  (right). 
Table  1:  MetalMapper  receiver  locations  with  respect  to  the  center  of  the  Z  transmitter  loop. 


Rx# 

X  [cm] 

Y  [cm] 

Z  [cm] 

0 

39 

39 

5 

1 

-26 

26 

5 

2 

13 

13 

5 

3 

0 

0 

5 

4 

-13 

-13 

5 

5 

26 

-26 

5 

6 

-39 

-39 

5 

The  MM  transmitters  are  modeled  as  infinitely  thin  rectangular  wires.  The  primary  magnetic 
induction  produced  at  any  observation  point  r  by  the  T-th  loop  is  determined  simply  from  the  Biot-Savart 
law, 


Br(r 

i= i 


-T,i 

Rl 


T  =  1,2,3, 


(162) 


where,  Rr  .  =  Ir -r'T  .1 ,  .is  the  location  of  the  z'-th  current  element,  and  A^  .  is  the  tangential  length 
vector  for  the  z'-th  subsection  of  the  loop.  In  what  follows,  and  unless  we  note  otherwise,  we  divide  each 
transmitter  coil  into  Nn  -  40  subsections  whenever  we  calculate  the  primary  magnetic  induction  using 
Eq.  (162). 
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Figure  8.  The  MetalMapper  geometry.  The  observation  point  r  is  defined  with  respect  to  the  global  Cartesian 
coordinate  system  XYZO;  r'  is  the  location  of  the  t-th  current  element  on  (in  this  case)  the  T=  3  transmitter,  which 

carries  a  current  /,  in  the  direction  p.  . 

3  -t  3,t 

The  MM  receiver  assembly  consists  of  seven  cube  sensors.  Each  of  these  measures  along  three 


orthogonal  directions  the  induced  voltages  that,  from  Faraday’s  law,  correspond  to  the  negative  of  the 
time  derivative  of  the  secondary  magnetic  flux  through  the  area  spanned  by  the  different  coils.  The 


induced  voltage  in  the  R-th  sensor  along  the  a- th  direction,  where  R  =  0,  ,6  and  a-z, y,x,  is 

computed  using 


(163) 


where  s“  is  the  area  of  the  relevant  coil  (all  of  which  are  10  cm  x  10  cm  squares  in  MetalMapper)  and 
n  is  the  unit  vector  perpendicular  to  it,  AsaR  and  r“K  are  respectively  the  i  -th  sub-area  and  vector 
location  point  on  s“  ,  B  is  the  magnetic  induction  (proportional  to  the  magnetic  field 

H  (r^))  produced  at  r"R  by  a  source  placed  at  r  .  Within  the  ONVMS  model,  Wir"K)  is  calculated 
using  equation  (68).  In  what  follows  we  always  divide  sa  into  N Rr  =  4  sub-areas. 
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Figure  9.  Response  of  an  81 -mm  mortar  illuminated  by  the  MM  Z-transmitter:  measured  (left),  ONVMS  prediction 
(center),  and  mismatch  between  modeled  and  actual  data  (right).  The  mortar  is  placed  35  cm  below  the  sensor  center 
and  oriented  45  degrees  nose  down.  The  data  are  plotted  in  logic  scale. 


To  validate  the  MetalMapper  versions  of  our  advanced  EMI  codes  we  conducted  comparisons 
between  actual  and  measured  data  for  different  targets.  Figure  9  through  Figure  1 1  compare  measured  and 
ONVMS-modeled  data  for  an  81 -mm  mortar  placed  35  cm  below  the  sensor  center,  oriented  45  degrees 
nose-down  and  illuminated  in  turn  by  the  Z,  Y,  and  X  transmitters.  We  use  three  responding  ONVMS 
sources  whose  locations  are  determined  with  the  combined  ONVMS -DE  algorithm.  The  inverted  location 
matches  the  actual  target  location  very  well.  The  model  is  seen  to  predict  target  EMI  responses  very 
accurately. 
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Figure  10.  Response  of  an  81 -mm  mortar  illuminated  by  the  MM  T-transmitter:  measured  (left),  ONVMS  prediction 
(center),  and  mismatch  between  modeled  and  actual  data  (right).  The  mortar  is  placed  35  cm  below  the  sensor  center 
and  oriented  45  degrees  nose  down.  The  data  are  plotted  in  logic  scale. 
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Figure  11.  Response  of  an  81 -mm  mortar  illuminated  by  the  MM  X-transmitter:  measured  (left),  ONVMS  prediction 
(center),  and  mismatch  between  modeled  and  actual  data  (right).  The  mortar  is  placed  35  cm  below  the  sensor  center 
and  oriented  45  degrees  nose  down.  The  data  are  plotted  in  log10  scale. 
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Figure  12:  Photo  of  the  TEMTADS  in  deployment  at  Blossom  Point  Test  Site  (left)  and  a  schematic  diagram  of  its 
Tx/Rx  sensors  (right). 


5.3  TEMTADS 

5.3.1  TEMTADS  modeling 

The  NRL  time-domain  EMI  sensor  array  TEMTADS  is  a  next-generation  system  designed  for 
subsurface  target  discrimination.  The  sensor  consists  of  25  transmit/receive  pairs,  each  composed  of  a  35- 
cm  square  transmitter  loop  surrounding  a  25-cm  square  receiver  loop,  arranged  in  a  rectangular  5x5  grid 
with  40-cm  neighbor-to-neighbor  separation  [83]  (Figure  12).  The  sensor  activates  the  transmitter  loops  in 
sequence,  one  at  a  time,  and  for  each  transmitter  all  receivers  receive,  measuring  the  complete  transient 
response  over  a  wide  dynamic  range  of  time  going  approximately  from  100  microseconds  (ps)  to  25 
milliseconds  (ms)  and  distributed  in  123  time  gates.  The  sensor  thus  provides  625  spatial  data  points  at 
each  location,  with  unprecedented  positional  accuracy. 

In  modeling  for  TEMTADS,  the  transmitter  loops  are  idealized  as  infinitesimally  thin 
35  cm  x  35  cm  square  loops.  The  primary  field  produced  at  any  observation  point  by  a  given  transmitter 
loop  is  determined  from  equation  (162).  We  use  NTx  =  20  for  TEMTADS  unless  we  note  otherwise.  The 
TEMATDS  measured  signal  is  modeled  using  equation  (163),  assuming  a  =  z  throughout  and  receiver 
sizes  of  25  cm  x  25  cm  and  dividing  each  receiver  into  NRx  =  9  sub-areas.  We  compare  actual  and 
ONVMS  modeled  data  for  a  105 -mm  projectile  in  Figure  13  and  find  very  good  agreement. 


Sky  Research,  Inc. 


59 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Figure  13.  Measured  (top  five  rows)  and  ONVMS-modeled  (bottom  five)  TEMTADS  data  for  a  105-mm  projectile 
at  the  25th  time  channel.  The  target  is  buried  at  a  depth  of  30  cm  and  oriented  horizontally  relative  to  the 
TEMATDS  system. 
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Figure  14:  The  APG  TOI. 


Table  2:  Inverted  location  and  orientation  for  TEMTADS  data 


Case  # 

Ground  truth  /estimated  for  a  37  mm  UXO 

1 

Xo  [m] 

Yo  [m] 

Zo  [m] 

Azimuth  [Degree] 

Dip  [Degree] 

2 

0.0/(0.03) 

0. 0/(0. 02) 

-0.35/(-0.39) 

0/(3) 

0/(5) 

3 

0.0/(0.013) 

0.0/(0.007) 

-0.34/(-0.369) 

0/(3) 

90/(88) 

4 

0.0/(0.001) 

0.0/(0.02) 

-0.38/(-0.41) 

0/(5) 

-90/(85) 

5 

0.0/(0.04) 

0.0/(0.05) 

-0.37/(-0.405) 

0/(5) 

45/(35) 

5.3.2  APG  test-site  classification 

To  demonstrate  the  classification  performance  of  the  advanced  EMI  models  we  conducted 
discrimination  studies  at  the  APG  test  site.  We  applied  a  combined  HAP/NSMS  approach  to  TEMTADS 
data  sets.  The  main  objective  of  the  study  was  to  discriminate  TOI  from  non-TOI  targets  and  further  to 
indicate  the  type  and  caliber  of  each  TOI.  The  TOI  at  APG  varied  in  size  from  25  mm  up  to  155  mm  and 
are  depicted  in  Figure  14. 

There  were  three  types  of  data  sets:  1)  Test  stand  data  set  collected  for  14  UXO  items  placed  in 
air  for  different  depths  and  orientations;  2)  Calibration  grid  data  sets  collected  over  the  same  targets  and 
over  some  clutter  items;  3)  Blind  grid  data  sets  collected  over  214  buried  items.  According  to  a 
preliminary  data  analysis  by  ESTCP,  soil  responses  were  insignificant  at  this  site,  and  they  were  thus 
subsequently  neglected.  The  test-stand  and  calibration  grid  data  sets  were  used  to  test  data  inversion  and 
discrimination  algorithms.  Object  depths  were  inverted  for  each  grid  using  the  HAP  method.  The  results 
for  a  37-mm  UXO  are  tabulated  in  Table  2.  Since,  TEMTADS  half  thickness  is  5cm,  the  inverted  depths 
were  in  very  good  agreement  (between  l=(-4+5)  and  2=(-3+5)  cm)  with  the  actual  depths  for  test-stand 
UXO  items. 
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Figure  15:  Comparison  between  the  inverted  and  actual  depth  for  all  65  APG  calibration  targets. 

We  also  used  the  HAP  method  to  invert  for  the  depths  of  all  65  calibration  targets.  The  results  are 
depicted  in  Figure  15.  The  inverted  depth  differed  by  up  to  15  cm  from  the  ground  truth,  a  difference  due 
to  the  fact  that  HAP  estimates  the  distance  from  the  sensor  center  to  the  target  center,  as  was  recorded  for 
test-stand  cases,  while  for  calibration  items  the  depths  were  measured  from  the  ground  surface.  The 
sensor  is  4"  (10.1  cm)  above  the  ground  and  the  transmitters  are  about  10  cm  thick,  and  therefore  the 
method  provides  reasonably  accurate  depth  estimates. 

Once  we  established  that  the  HAP  method  estimates  depths  accurately  for  test-stand  and 
calibration  items  we  proceeded  to  estimate  the  total  NSMS  for  all  items  and  used  it  for  discrimination. 
Figure  16  shows  the  inverted  total  NSMS  as  a  function  of  time  from  test-stand  TEMTADS  data  sets  with 
the  105-mm  projectile  and  the  81 -mortar  as  targets.  Each  set  of  test-stand  measurements  comprised  six 
different  depths  and  target  orientations.  The  total  NSMS  is  seen  to  be  unique  for  all  cases  and,  for  both 
test-stand  and  calibration  data.  We  then  determined  the  best  NSMS  classification  features.  We  fit  the  total 

M__  NSMS  curves  with  the  Pasion-Oldenburg  expression  Mzz(t )  =  kt  ,  where  t  is  time,  and  k  ,  [>  and 

y  are  the  fitting  parameters  for  each  anomaly.  We  studied  different  combinations  of  In  k  ,  ft ,  and  y  using 
test-stand  data.  The  results  for  f5  vs.  Ink  appear  in  Figure  17,  and  those  for  y  vs.  Ink  and  y  vs.  ft 
appear  in  Figure  16.  We  see  that  the  best  classification  performance  is  achieved  using  Ink  and  f3 . 
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Figure  16:  Inverted  total  NSMS  for  APG  test-stand  105  mm  projectile  and  81  mm  mortar. 
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Figure  17:  Scatter  plot  of  inverted  /?  vs.  In  k  classification  features  for  APG  test-stand  TOI. 
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Figure  18:  Scatter  plot  of  inverted  y  vs  In  k  (left)  and  (i  (right)  parameters  for  APG  test-stand  TOI. 
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Figure  19:  Scatter  plot  of  inverted  /?  vs.  In  k  classification  features  for  all  214  APG  blind-test  anomalies. 


Finally,  the  described  data  inversion  and  classification  schemes  were  applied  to  the  214  blind 
grid-data  cells.  These  were  also  first  inverted  to  determine  the  total  NSMS,  from  which  time -decay- 
history  curves  were  synthesized,  discrimination  features  were  extracted,  and  classification  was  performed 
vis-a-vis  test-stand  UXO  items.  A  scatter  plot  of  inverted  and  classified  In  k  and  /?  features  for  all  214 
APG  test  anomalies  is  shown  in  Figure  19.  The  result  illustrates  that  the  inverted  features  for  60-mm,  81- 
mm,  and  105-mm  TOI  are  clustered  tightly,  while  those  for  37-mm  and  25-mm  TOI  s  are  scattered  and 
mixed  with  those  of  clutter  items.  This  complicates  classification. 

To  overcome  this  problem,  in  addition  classification/clustering  approach,  the  entire  time  decay 
history  of  the  total  NSMS  were  also  examined  and  compared  to  the  total  NSMS  of  the  test-stand  TOI 
case-by-case  as  a  check  on  the  classification.  The  comparisons  are  summarized  in  Figure  20  and  Figure 
2 1 .  For  all  APG  test  anomalies  a  ranked  list  was  created  in  which  the  anomalies  were  ranked  as  clutter  or 
TOI  and  TOI  were  further  ranked  by  caliber  and  type.  This  list  was  submitted  to  the  Institute  for  Defense 
Analyses  (IDA)  for  independent  scoring.  The  scores  showed  that  the  advanced  model  was  able  to  identify 
all  UXO  as  TOI  and  classified  all  UXO  correctly  by  type  and  caliber.  The  false-positive  rate  was  5%. 
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Figure  20:  Comparison  between  library  (green  lines)  and  inverted  (red  and  blue  lines)  blind-test  total  NSMS  for 
105-mm  projectiles,  81-mm  munitions,  and  60-mm  mortars. 
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Figure  21:  Comparisons  between  library  (green  lines)  and  inverted  (red  and  blue  lines)  blind-test  total  NSMS  for 
37-mm  and25-mm  mortars. 
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Figure  22:  Schematic  diagram  of  the  BUD  system. 


+X 


5.4  BUD 

The  Berkeley  UXO  discriminator  (BUD)  is  an  advanced  standalone  time-domain  system 
developed  at  the  University  of  California  to  detect  and  discriminate  UXO  in  the  20-mm  to  155-mm  size 
range,  and  consists  of  three  orthogonal  coil  transmitters.  The  horizontal  Z-coils  are  vertically  separated  by 
26"  and  have  a  39"  x  39"  footprint.  The  Y-  and  X-vertical  coils  are  mounted  on  the  diagonals  between  the 
Z-coils  (see  Figure  22):  the  X-coils  are  45.5"  x  23.5"  while  the  F-coils  are  45.5"  x  22.5"  in  size,  and  both 
are  separated  by  6".  The  BUD  illuminates  targets  in  three  independent  directions,  which  induce  eddy 
currents  in  all  three  modes.  BUD  has  eight  pairs  of  differenced  receiver  coils  placed  horizontally  along 
the  two  diagonals  of  the  upper  and  lower  planes  of  the  Z-transmitter  loops.  The  pairs  are  located  on 
symmetry  lines  trough  the  center  and  are  wired  in  opposition  so  as  to  cancel  the  primary  magnetic  field 
during  transmission.  Figure  22  shows  the  BUD  system  in  operation. 

The  BUD  transmitter  loops  were  modeled  as  idealized  infinitely  thin  square  loops.  The  primary 
fields  produced  at  any  observation  point  by  the  transmitters  are  determined  using  a  suitable  modification 
of  equation.  (162),  again  with  N  =  40 .  The  BUD  measured  signals  are  modeled  using  equation  (163)  as 


1 'KX 

v*  =  -I 


dB, -(rf.*  ~ro) 
dt 


Nrx  (r»r  — r  ^ 

*»,.,+£  ,QV  °: V. 


dt 


As; 


=  A^z 


(164) 


where  r  and  r!  are  the  locations  of  the  Rx  and  Rx'  receivers,  given  in  Table  3.  For  the  case  of  BUD 
we  divide  the  receivers  into  NB  =  9  sub-areas. 

Rx 
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Figure  23:  The  BUD  system  in  operation. 


Table  3:  BUD  receiver  locations  with  respect  to  the  origin. 


Rx# 

X  [cm] 

Y  [cm] 

Z  [cm] 

Rx'# 

X'[cm] 

l"[cm] 

Z'[cm] 

1 

35.48 

35.48 

0 

r 

-35.48 

-35.48 

66 

2 

-35.48 

35.48 

0 

T 

35.48 

-35.48 

66 

3 

-35.48 

-35.48 

0 

3’ 

35.48 

-35.48 

66 

4 

35.48 

-35.48 

0 

4’ 

-35.48 

35.48 

66 

5 

19.29 

19.29 

0 

5’ 

-19.29 

-19.29 

66 

6 

-19.29 

19.29 

0 

6’ 

19.29 

-19.29 

66 

7 

-19.29 

-19.29 

0 

T 

19.29 

19.29 

66 

8 

19.29 

-19.29 

0 

8’ 

-19.29 

19.29 

66 

All  data  presented  here  were  collected  by  personnel  from  the  Berkeley  UXO  team  at  Yuma 
Proving  Ground  in  Arizona  over  objects  at  different  orientations  and  depths.  The  response  of  each  object 
was  represented  with  only  five  NSMS.  Figure  24,  Figure  25,  and  Figure  26  show  comparisons  between 
modeled  and  actual  data  for  all  transmitters  and  receivers  and  for  all  time  channels.  The  results  clearly 
show  that  the  NSMS  very  well  predicts  the  EMI  response  of  a  M-75  mm  UXO.  Total  NSMS  amplitudes 
were  determined  for  three  samples  each  of  M-75,  60-mm,  and  M-37  UXO  and  are  depicted  in  Figure  27. 
The  result  demonstrates  that  the  NSMS  is  applicable  to  the  BUD  system  and  is  a  good  discriminator. 
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Figure  24:  Comparisons  between  actual  and  predicted  data  for  an  M75  UXO  illuminated  by  the  BUD  Z  transmitter. 
Solid  lines  are  actual  data,  circles  stand  for  NSMS  predictions. 


Figure  25:  Comparisons  between  actual  and  predicted  data  for  an  M75  UXO  illuminated  by  the  BUD  X  transmitter. 
Solid  lines  are  actual  data,  circles  stand  for  NSMS  predictions. 
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Figure  26:  Comparisons  between  actual  and  predicted  data  for  an  M75  UXO  illuminated  by  the  BUD  Y  transmitter. 
Solid  lines  are  actual  data,  circles  stand  for  NSMS  predictions. 
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Figure  27:  Recovered  total  NSMS  from  calibration  BUD  measurements  for  M-75  (blue),  37-mm  (green),  and 
M-60  (red)  UXO. 
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Figure  28.  Photo  and  schematic  diagram  of  the  MPV  sensor. 


5.5  MPV 

The  MPV  sensor,  developed  by  G&G  Sciences,  Inc.,  consists  of  two  transmitter  loops  and  five 
triaxial  receiver  cubes.  The  receivers  are  located  as  follows:  Cube  #0  above  center  (z  -  30.6  cm);  Cube  #1 
at  the  origin;  Cube  #2  left  of  center  (x  =  -39.6  cm);  Cube  #3  forward  of  center  (y  =  39.6  cm);  and  Cube  #4 
right  of  center  ( x  =  39.6  cm).  These  receivers  accurately  measure  the  complete  transient  response  over  a 
wide  dynamic  range  of  time  going  from  100  ps  to  25  ms.  In  numerical  models  we  assume  that  the 
transmitter  loops  are  idealized  as  infinitely  thin  circular  loops  with  37.5  cm  radii,  and  separated  by  12  cm. 
The  complete  primary  field  produced  at  any  observation  point  by  the  transmitter  loop  is  determined  from 
equation  (162)  as 


B(r)  =  ^XZ/A*,,<3XR,'‘ 

4/rtrtr  r. 


(165) 


where,  for  the  Mh  transmitter  loop,  t  —  1,2,  R  =  r-r'.  ,  r'ti  is  the  location  of  the  /-th  current  element 

on  t-transmitter,  and  Af  ( .  is  the  tangential  length  vector  for  the  i- th  subsection.  We  use  N  =  20  unless 
we  note  otherwise.  The  MPV  measured  signal  is  modeled  using  equation  (163)  with  each  loop  having 
area  10  x  10  cm2  and  divided  into  4  sub-areas. 
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Figure  29:  Multi-object  MPV  data  collection  setup  (right).  The  red  circle  corresponds  to  the  MPV  head,  which  was 
placed  stationary;  the  targets  were  moved  along  the  blue  line.  The  center  of  the  first  target  (the  81 -mm)  was  placed 
at  the  blue  points,  and  the  distance  between  the  first  and  second  targets  was  kept  fixed. 

To  illustrate  the  applicability  of  the  ONVMS  for  MPV  data  we  conducted  studies  in  multi-target 
inversion  and  discrimination.  The  measurements  reported  here  were  conducted  at  the  SKY  Research 
office  in  Hanover,  New  Hampshire.  The  sensor  was  placed  stationary,  and  data  were  collected  for  two 
objects  with  different  separations  and  orientations  placed  on  5  x  5  grid  points.  The  separation  between  the 
grids  points  was  20  cm.  The  targets  were  an  8 1  -mm  munition  and  a  40-mm  round.  The  data  were  inverted 
using  the  simple  dipole  model  with  DE  and  the  ortho-normalized  volume  magnetic  source  model 
(ONVMS).  The  number  was  assumed  given  in  the  simple  dipole  model,  while  in  the  ONVSMS  four 
arbitrarily  distributed  interacting  dipoles  were  used.  The  dipoles’  positions  were  determined  using  DE. 
The  inverted  polarizability  tensor  principal  elements  for  the  projectiles  are  depicted  in  Figure  30  for  three 
different  target-to-target  separation  vectors:  (-25,  0,  0)  cm  (blue),  (-40,  0,  0)  cm  (red),  and  (-25,  0, 
25)  (green).  The  single-dipole/DE  algorithm  accurately  inverts  the  polarizability  elements  for  the  shallow 
81 -mm  projectile  but  fails  to  identify  the  40  mm  projectile  when  the  distance  between  the  two  is  25  cm 
(blue)  and  when  the  40-mm  is  placed  deeper  (green).  When  the  distance  between  the  targets  increases  and 
they  both  have  the  same  depth  the  algorithm  identifies  the  40-mm  projectile  correctly.  The  same  data  sets 
were  inverted  using  the  combined  ONVMS -DE  technique.  The  inverted  locations  showed  the  ONVMS 
dipoles  grouped  around  the  locations  of  the  projectiles,  and  for  discrimination  we  summed  the  ONVMS 
amplitudes  for  each  group.  The  results  for  the  two  targets,  which  appear  in  Figure  31,  show  that  the 
inverted  ONVMS  is  consistent  for  all  cases  and  both  munitions.  The  ONVMS  technique  is  seen  to  be  a 
robust  algorithm  for  discriminating  not  only  single  well-isolated  targets  but  also  multi-target  scenarios. 
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Figure  30:  Inverted  polarizability  principal  elements  for  two  targets  in  three  different  setups;  results  for  the  81 -mm 
projectile  at  left  and  for  the  40-mm  munition  at  right.  In  all  three  cases  the  targets  were  horizontal,  and  the  vertical 
distance  between  the  MPV  center  and  the  81 -mm  was  40  cm.  The  center  to  the  center  coordinate  differences 
between  the  81 -mm  and  40-mm  projectiles  are  (-25,  0,  0)  cm,  (-40,  0,  0)  cm,  and  (-25,  0,  25)  cm. 


Figure  31:  Inverted  total  ONMS  for  81  mm  (left)  and  40  mm  (right)  projectiles  for  three  different  cases. 

We  have  just  compared  the  single -dipole  and  ONVMS  model  for  UXO  discrimination.  (We  do 
note  that  in  all  cases  we  used  DE  to  perform  the  crucial  task  of  determining  object  locations).  The  dipole 
model  is  sufficient  for  inversion  when  the  different  targets  are  well  separated  but  breaks  down  when  they 
are  placed  close  to  each  other  or  when  the  EMI  response  from  one  item  dominates.  In  contrast,  the 
physically  complete  model  is  able  to  predict  target  EMI  responses  accurately  for  these  situations,  making 
the  ONVMS  method  our  preferred  tool  for  the  live-site  UXO  classification  studies  we  present  next. 
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6  ESTCP  live-site  classification  studies  using  advanced  models 

6.1  Introduction 

The  Environmental  Security  Technology  Certification  Program  (ESTCP)  recently  launched  a 
series  of  live-site  UXO  classification  blind  tests  at  increasingly  challenging  and  complex  sites  [84-86] 
aiming  to  demonstrate  the  performance  of  advanced  EMI  detection  technologies  and  UXO  discrimination 
and  classification  algorithms.  The  first  test  was  conducted  in  2007  at  the  UXO  live  site  at  the  former 
Camp  Sibert  in  Alabama  using  first-generation  EMI  sensors  (the  commercially  available  EM61-MK2  and 
EM-63,  both  developed  by  Geonics  Ltd.).  The  Sibert  test  was  relatively  simple:  one  had  to  discriminate 
well-isolated  large  intact  4.2"  mortars  from  smaller  range  scrap,  shrapnel,  and  cultural  debris.  The  second 
ESTCP  discrimination  study  to  demonstrate  the  applicability  of  EMI  classification  technologies  was  set 
up  in  2009  at  the  live  UXO  site  in  San  Luis  Obispo  (SLO)  in  California  and  featured  a  more  challenging 
topography  and  a  wider  mix  of  TOI  [84-85].  Magnetometers  and  first-generation  EMI  sensors  were 
deployed  on  the  site  and  used  in  survey  mode.  Two  advanced  EMI  sensing  systems — the  Berkeley  UXO 
Discriminator  (BUD)  of  Section  5.4  and  the  Naval  Research  Laboratory’s  TEMTADS  EMI  array, 
presented  in  Section  5.3 — were  used  to  perform  cued  interrogation  of  the  anomalies  detected.  A  third 
advanced  system,  the  Geometries  MetalMapper  of  Section  5.2,  was  used  in  both  survey  and  cued  modes 
for  identifying  and  classifying  anomalies.  Among  the  munitions  buried  at  SLO  were  60-mm  and  81 -mm 
projectiles,  4.2"  mortars,  and  2.36"  rockets;  three  additional  munition  types  were  discovered  during  the 
course  of  the  demonstration.  The  third  site  chosen  was  the  former  Camp  Butner  in  North  Carolina.  That 
demonstration  was  designed  to  investigate  evolving  classification  methodologies  at  a  site  contaminated 
with  37-mm  projectiles,  adding  yet  another  layer  of  complexity  into  the  process  [87-89].  In  this  chapter 
we  describe  the  work  we  performed  when  we  participated  in  those  studies  and  summarize  the  results  we 
obtained. 

6.2  Camp  Sibert 

In  2006,  researchers  affiliated  with  Sky  Research,  Inc.  collected  data  at  Camp  Sibert  using  the 
EM-63,  a  cart -based  step-off  time-domain  EMI  sensor  produced  by  Geonics  Ltd.  [90].  The  targets  buried 
in  216  cells — some  of  which  were  empty — included  unexploded  4.2"  mortar  shells,  mortar  explosion 
byproducts  like  base  plates  and  partial  mortars  (i.e.,  stretched-out  half-shells),  smaller  shrapnel,  and 
unrelated  metallic  clutter;  some  examples  appear-  in  Figure  32.  The  different  items  were  distributed  as 
shown  in  Figure  32(d). 
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(a)  4.2"  mortar  shell 


(b)  Base  plate 


Type 

Training 

Testing 

Total 

UXO 

38 

34 

72 

Partial 

12 

23 

35 

Base 

5 

40 

45 

Scrap 

6 

25 

31 

Clutter 

4 

22 

26 

Empty 

1 

6 

7 

Total 

66 

150 

216 

(c)  Half-shell 


(d)  Cell  contents 


Figure  32:  Camp  Sibert  anomalies:  4.2  inch,  base  plates  and  partial  mortars. 

We  analyzed  the  Sibert  data  using  HAP  and  NSMS.  By  combining  those  two  techniques  we  made 
sure  our  method  of  analysis  [82,  91]  avoided  the  tendency  of  inversion  algorithms  to  linger  in  local 
minima.  We  performed  the  localization  step  independently  at  the  outset  and  then  used  its  results  to  help  in 
the  characterization,  allowing  for  fast  and  accurate  determination  of  the  total  NSMS  for  each  target.  We 
classified  these  NSMS  values  using  a  heuristic  pattern-matching  method  (Section  6.2.1),  an  open-source 
implementation  [92]  of  SVM  (Section  6.2.3),  and  mixed  modeling  (Section  6.2.4).  The  SVM-based 
classification  improved  upon  template-matching  [93-94]  in  that  it  required  less  human  intervention  and 
was  thus  faster  to  run  and  easier  to  adapt  to  other  sets  of  observations.  On  the  other  hand,  the  semi- 
supervised  Gaussian  mixture  model  provided  a  classification  performance  exceeding  that  of  SVM,  which 
made  it  our  preferred  statistical  classification  procedure  for  use  in  all  subsequent  classification  tasks. 

6.2.1  Target  location  and  characterization;  preliminary  pattern-matching  classification 

We  started  the  procedure  by  applying  HAP  to  determine  the  target  location  for  each  cell.  Figure 
33  compares  actual  and  inverted  data  at  the  first  and  20th  time  channels  (top  and  bottom  rows 
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respectively)  for  one  cell.  To  find  the  target  we  take  a  fictitious  5  m  x  5  m  flat  square  surface  concentric 
with  the  plot  and  located  30  cm  below  the  sensor  (i.e.,  at  ground  level)  and  divide  it  into  11x11  patches, 
each  of  which  is  assumed  to  contain  a  magnetic -charge  distribution  of  uniform  density.  We  take  the 
measured  field  data  (seen  on  the  left  column  of  Figure  33)  and  use  Eq.  (109)  to  determine  q  ,  which  in 
turn  allows  us  to  determine  y/(r)  using  Eq.  (Ill)  and  construct  the  matrices  of  Eq.  (108)  to  find  the 
location.  We  do  this  separately  for  every  time  channel  and  get  consistent  location  estimates  from  gate  to 
gate,  which  lends  credence  to  their  precision.  The  depths  thus  determined  are  also  acceptably  close  to  the 
ground  truth. 

After  finding  the  locations  we  run  a  fully  three-dimensional  orientation-free  NSMS  code  to 
determine  the  time -dependent  total  NSMS  amplitude  for  all  cells.  To  compute  <2(0  we  surround  the 
target  with  a  prolate  spheroid  of  semiminor  axis  a  —  5  cm  and  elongation  e  =  b  /  a  -  4 .  This  spheroid  is 
divided  into  seven  azimuthal  belts,  each  of  which  is  assumed  to  contain  a  radial-magnetic-dipole 
distribution  of  constant  density.  The  spheroid  is  placed  at  the  location  estimated  by  the  HAP  method  and 
the  orientation  given  by  the  dipole  moment  m  obtained  from  Eqs.  (106)  and  (99)-(101).  With  all  the 
pieces  in  place,  we  apply  Eq.  (25)  to  find  Q  and  Eq.  (26)  to  extract  Q(t)  for  the  target.  The  inverted  total 

NSMS  for  all  anomalies,  and  for  4.2"  mortars,  base  plates,  and  partial  mortars  are  depicted  in  Figure  34. 

It  is  evident  that  there  are  distinguishable  differences  between  the  total  NSMS  for  the  4.2" 
mortars,  the  base  plates,  and  the  partial  mortars.  Particularly  at  late  times,  each  target  has  different  natural 
decay  characteristics  that  depend  on  its  geometry  and  material  properties.  It  is  also  important  to  notice 
that  the  total  NSMS  for  the  4.2”  mortars  is  very  well  grouped.  To  further  simplify  the  classification  task 
we  used  the  Pasion-Oldenburg  law  to  fit  the  time-dependent  NSMS  curves,  obtaining  as  a  result  the 
amplitudes  ( k  ),  the  power-law  exponents  ( ft ),  and  the  exponential-decay  inverse  time  constants  ( ;/  ),  all 
of  which  we  tested  as  classification  features.  We  obtained  the  parameters  by  direct  nonlinear  least-squares 
fit  of  (57)  and  by  linear  (pseudo)inversion  of  its  logarithm  (166);  both  procedures  gave  consistent  results. 
In  general  we  obtain  good  fits  to  the  measured  fields  [94];  Figure  33  shows  that  the  discrepancy  between 
the  actual  data  and  the  model  prediction  runs  only  to  a  few  percent. 
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Figure  33:  Camp  Sibert  EM-63  near  field  distributions:  Left  and  middle  columns:  actual  and  modeled  data 
respectively.  Right  column:  misfits. 

After  investigating  different  combinations  of  these  feature-space  parameters  we  found  that  k  in 
conjunction  with  the  ratio  Q(tl5)  /  Q(t{)  which  involves  a  fixed  superposition  of  ft  and  y  ,  worked  best: 

the  left  panel  of  Figure  35  depicts  this  winning  combination  for  all  items  and  clearly  shows  the  tight 
clustering  and  generous  cluster-to-cluster  separation  that  generally  lead  to  reliable  classification.  (The 
15th  time  channel,  centered  at  about  2.7  ms,  was  chosen  because  it  takes  place  late  enough  to  show  the 
behavior  described  above  but  early  enough  that  all  targets  still  have  an  acceptable  signal-to-noise  ratio; 
nearby  time  channels  produce  similar  results.)  When  we  received  the  ground  truth  for  all  targets  we 
proceeded  to  construct  the  ROC  curve  that  appears  in  the  right  panel  of  Figure  35.  We  see  that  only  one 
excavation  out  of  130  anomalies  is  necessary  before  all  UXO  are  indentified  correctly. 

We  obtain  similar  results  using  the  SVM  algorithm. 


Sky  Research,  Inc. 


77 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


o 


o 


Figure  34:  Inverted  total  NSMS  for  all  anomalies:  4.2"  mortars,  base  plates,  and  partial  mortars. 
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Figure  35:  Left:  Classification  features.  Right:  ROC  curve  of  NSMS  performance. 
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(a)  Shell  at  No.  7  (b)  Twisted  chain  (c)  Crumpled  debris 

Figure  36:  a)  Unexploded  shell  from  Cell  No.  7  and  (b,c)  the  two  false  alarms  obtained  by  the  SVM  classifier  using 
k  and  <2(f15 )  /  Q(t] )  as  discriminators. 


6.2.2  SVM  classification 

We  use  a  Gaussian  RBF  kernel  for  the  SVM  analysis.  The  kernel  width  turns  out  not  to  have 
much  influence  on  the  outcome;  we  usually  set  it  so  that  a  unit  in  a  typical  x-  or  y-axis  in  a  log  plot  (for 
example,  Figure  37)  comprises  100A  Gaussian  widths,  where  A  is  the  dimensionality  of  the  feature  space. 
To  find  the  capacity  C  we  train  the  SVM  with  a  subset  of  the  training  data  and  a  given  C,  scramble  the 
training  set,  and  use  a  new  subset  of  the  data  for  testing.  We  then  vary  C,  setting  it  to  a  high  value  initially 
and  then  lowering  it,  and  keep  the  lowest  capacity  with  which  the  machine  identifies  all  dangerous  items 
in  the  test.  The  procedure  is  rather  ad  hoc  but  effective  for  the  data  at  hand,  given  the  small  sample  sizes, 
the  low  dimensionality  of  the  feature  spaces,  and  the  speed  of  the  SVM  implementation.  A  more 
systematic  search  for  C  and  y  using  five-fold  cross-validation  [92]  recommends  slightly  higher  capacities 
that  result  in  identical  predictions. 

For  R  and  k  as  features  we  find  the  best  SVM  performance  using  C  -  10.  The  results  are 
displayed  (for  testing  data  only)  in  Table  4  and  shown  pictorially  (for  both  training  and  testing)  in  Figure 
37.  The  matrix  element  c..  in  the  table  denotes  an  item  of  category  i  that  was  identified  by  the  SVM  as 

belonging  to  category  y;  in  other  words,  the  rows  of  this  contingency  table  correspond  to  the  ground  truth 
and  the  columns  to  predictions.  The  small  markers  in  the  plot  show  the  ground  truth  (hollow  for  training 
data  and  filled  for  the  tests),  while  the  large  markers  point  out  the  items  for  which  the  SVM  makes  wrong 
predictions.  For  example,  a  small  yellow  upright  triangle  surrounded  by  a  large  cyan  square  is  a  piece  of 
scrap  (clutter  unrelated  to  UXO)  incorrectly  identified  as  a  base  plate.  The  UXO,  with  their  high  initial 
amplitudes  and  slow  decay,  are  clustered  at  the  top  right  corner.  We  see  that  there  are  only  two  false 
alarms  (i.e.,  objects  identified  with  UXO  that  were  in  fact  something  else)  and  that  all  potentially 
dangerous  items  have  been  identified  correctly. 


Sky  Research,  Inc. 


79 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


log  R 


Figure  37:  Result  of  the  SVM  classification  for  the  Camp  Sibert  anomalies  using  the  logarithms  of  k 
and  R  =  Q(tl5) /  Q(t{) .  The  SVM  has  been  trained  with  capacity  C  =  10  and  kernel  width  a  =  1/200.  The  small 

markers  denote  the  ground  truth  for  both  training  (hollow)  and  testing  (solid)  cells.  The  larger  markers  highlight  the 
cases  where  there  is  disagreement  between  the  ground  truth  and  the  SVM  prediction. 

Table  4:  SVM  classification  of  Camp  Sibert  anomalies  using  k  and  R  with  C  =  10 


SVM  prediction 


k, 

R-.  C  =  10 

UXO 

Partial 

Base 

Clutter 

Scrap 

Empty 

UXO 

34 

0 

0 

0 

0 

0 

3 

S— • 

Partial 

0 

22 

0 

1 

0 

0 

T3 

Base 

0 

0 

39 

1 

0 

0 

3 

Clutter 

0 

0 

4 

19 

0 

2 

o 

b 

C1 

Scrap 

2 

0 

3 

4 

13 

0 

Empty 

0 

1 

1 

1 

2 

1 

The  false  alarms,  two  pieces  of  non-UXO  clutter,  appear  in  Figure  36(b)  and  Figure  36(c).  They 
are  seen  to  be  similar  to  the  4.2"  mortars  in  size  and  metal  content  (cf.  Figure  36  (a)),  which  makes  their 
k  and  R  values  lie  closer  to  the  tight  UXO  cluster  than  to  any  other  anomaly.  Here  we  note  that,  as  can 
be  seen  in  Figure  32(d),  the  training  data  provided  by  the  examiners  was  somewhat  biased  toward  UXO, 
while  clutter  and  scrap  samples  were  underrepresented  (this  was  not  the  case  with  the  testing  data  and 
should  not  be  expected  in  future  tests).  If  we  switch  training  and  testing  data  in  the  SVM  analysis  we  can 
achieve  perfect  discrimination  without  varying  the  capacity — though  in  this  case  we  have  more  training 
data  than  tests.  This  highlights  the  importance  of  having  a  diverse  collection  of  representative  samples  to 
use  during  the  training  stage. 
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Table  5:  SVM  classification  of  Camp  Sibert  anomalies  using  y  and  k  with  C  =  9 


SVM  prediction 


UXO 

Partial 

Base 

Clutter 

Scrap 

Empty 

UXO 

34 

0 

0 

0 

0 

0 

Partial 

5 

17 

0 

1 

0 

0 

Base 

0 

0 

39 

0 

1 

0 

Clutter 

0 

0 

4 

15 

5 

1 

Scrap 

2 

1 

3 

5 

11 

0 

Empty 

1 

1 

2 

2 

0 

0 

We  can  repeat  the  analysis  using  other  two-dimensional  combinations  of  the  Pasion-Oldenburg 
parameters.  Combining  k  and  y  yields  results  similar  to  those  of  k  and  R ,  as  Figure  38  and  Table  5  show. 
Figure  39  and  Table  6  show  the  classification  resulting  from  the  use  of  (>  and  y  as  discriminators.  The 
table  shows  that  we  can  obtain  reasonable  discrimination,  with  all  the  UXO  once  again  correctly 
identified,  but  the  increased  number  of  false  alarms  and  the  very  high  capacity  needed  (four  orders  of 
magnitude  larger  than  the  previous  ones)  indicate  that  this  combination  of  parameters  may  not  be  optimal 
and  that  this  machine  is  prone  to  overfitting.  A  glance  at  the  figure  shows  the  clustering  is  much  less 
clear-cut  than  in  the  previous  cases,  partly  because  the  range  of  [>  is  rather  small.  In  fact,  combining  k  and 
fi  greatly  reduces  the  performance,  since  the  small  /l- range  and  the  close  similarity  in  k  of  the  UXO  and 
the  partial  mortars  cause  an  overlap  between  the  two  categories  that  cannot  be  disentangled. 

It  is  helpful  and  straightforward  to  increase  the  dimensionality  of  the  feature  space.  Figure  40 
shows  the  discrimination  obtained  by  running  the  SVM  using  all  three  Pasion-Oldenburg  features.  The 
capacity  C  =  9  here,  and  increasing  it  changes  the  results  only  slightly.  The  number  of  false  alarms 
increases:  we  get  the  same  two  pieces  of  scrap  from  before,  and  now  a  few  of  the  partial  mortars  are 
identified  as  UXO  by  the  algorithm,  due  in  part  to  the  small  range  of  and  in  part  to  the  large  gap 
between  the  UXO  and  the  other  anomalies,  clearly  visible  in  the  figure,  which  again  calls  out  for  more 
and  more -diverse  training  information. 

Finally,  it  is  possible  to  dispense  with  the  Pasion-Oldenburg  model  altogether  and  run  an  SVM 
using  the  “raw”  Q{t )  as  input.  The  feature  space  has  dimensionality  A  =  25.  We  scale  the  values  by 

<2(?j)  and  take  the  logarithm.  We  find  C  =  20  to  be  the  optimal  value.  Table  4  shows  the  results.  The 

performance  is  slightly  inferior  to  that  of  R  vs.  k:  the  usual  two  false  alarms  are  there,  along  with  a  few 
new  ones.  All  the  UXO  are  identified  correctly.  We  can  also  use  the  logarithm  of  Q  without  any  scaling 
(though  the  SVM  internally  rescales  the  feature  space  to  [0,1]A).  A  capacity  C  =  1  suffices  here.  The 
results  appear  on  Table  5.  All  dangerous  items  are  once  more  identified  as  such. 
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Figure  38:  Result  of  the  SVM  classification  for  the  Camp  Sibert  anomalies  using  the  logarithms  of  the  Pasion- 
Oldenburg  parameters  k  and  y.  The  SVM  here  has  a  capacity  C  =  9.  The  small  markers  denote  the  ground  truth  for 
both  training  (hollow)  and  testing  (solid)  cells.  The  larger  markers  show  the  wrong  SVM  predictions. 
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Figure  39:  Result  of  the  SVM  classification  for  the  Camp  Sibert  Anomalies  using  the  logarithms  of  the  Pasion- 
Oldenburg  parameters  j3  and  y.  The  SVM  capacity  C  =  105.  The  small  markers  denote  the  ground  truth  for  both 
training  (hollow)  and  testing  (solid)  cells.  The  larger  markers  highlight  the  wrong  predictions  made  by  the  SVM. 
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Figure  40:  SVM  classification  of  the  Camp  Sibert  Anomalies  using  the  logarithms  of  k,  b,  and  g.  The  SVM  has 
C  =  9.  The  small  markers  denote  the  ground  truth  for  both  training  and  testing  cells.  The  larger  markers  highlight  the 
cases  where  there  is  disagreement  between  the  ground  truth  and  the  SVM  prediction. 


Table  6:  SVM  classification  of  Camp  Sibert  anomalies  using  P  and  y  with  C  =  105 


SVM  prediction 


UXO 

Partial 

Base 

Clutter 

Scrap 

Empty 

UXO 

34 

0 

0 

0 

0 

0 

Partial 

0 

14 

5 

2 

1 

1 

Base 

3 

0 

37 

0 

0 

0 

Clutter 

0 

1 

5 

14 

1 

4 

Scrap 

3 

1 

1 

3 

13 

1 

Empty 

2 

0 

0 

3 

1 

0 

Table  7:  SVM  classification  of  Camp  Sibert  anomalies  using  the  complete  NSMS  time  decay 


SVM  prediction 


UXO 

Partial 

Base 

Clutter 

Scrap 

Empty' 

g 

UXO 

34 

0 

0 

0 

0 

0 

9 

Partial 

0 

15 

0 

*T 

/ 

1 

0 

-a 

Base 

3 

0 

34 

3 

0 

0 

c 

Clutter 

0 

2 

3 

14 

4 

2 

o 

i-t 

O 

Scrap 

3 

1 

3 

3 

12 

0 

Empty' 

2 

2 

1 

0 

1 

0 
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6.2.3  SVM  analysis  of  Camp  Sibert  data:  summary 

In  this  section  we  applied  the  NSMS  model  to  EM-63  Camp  Sibert  discrimination  data.  First  the 
locations  of  the  objects  were  inverted  for  by  the  fast  and  accurate  dipole-inspired  HAP  method. 
Subsequently  each  anomaly  was  characterized  at  each  time  channel  through  its  total  NSMS  strength. 
Discrete  intrinsic  features  were  selected  and  extracted  for  each  object  using  the  Pasion-Oldenburg  decay 
law  and  then  used  as  input  for  a  support  vector  machine  that  classified  the  items. 

Our  study  reveals  that  the  ratio  of  an  object’s  late  response  to  its  early  response  can  be  used  as  a 
robust  discriminator  when  combined  with  the  Pasion-Oldenburg  amplitude  k  .  Other  mixtures  of  these 
parameters  also  result  in  good  classifiers.  Moreover,  we  can  use  Q  directly,  completely  obviating  the 
need  for  the  Pasion-Oldenburg  fit.  In  each  case  the  classifier  runs  by  itself  and  does  not  require  any 
human  intervention.  The  SVM  can  be  trained  very  quickly,  even  when  the  feature  space  has  more  than  20 
dimensions,  and  it  is  a  simple  matter  to  add  more  training  data  on-the-fly.  It  is  also  possible  to  use  already 
processed  data  to  classify  examples  as  yet  unseen. 

We  should  stress  that  none  of  our  classifications  yielded  false  negatives:  all  UXO  were  identified 
correctly  in  every  instance.  (This  is  due  in  part  to  the  clean,  UXO-intensive  training  data  provided  by  the 
examiners  and  may  change  under  different  conditions.)  The  number  of  false  alarms  (false  positives) 
varies  with  the  classification  features,  but  is  in  general  low  and  can  be  as  low  as  2  out  of  36  reported 
positives.  Figure  37,  Figure  38  and  Figure  39  show,  among  others,  how  these  false  alarms  occur:  Some  of 
the  clutter  items  have  a  response  that  closely  resembles  that  of  UXO.  While  this  will  inevitably  arise,  it 
may  still  be  possible  to  make  the  SVM  more  effective — and  perhaps  approach  100%  accuracy — by 
including  some  of  these  refractory  cases  during  the  training.  That  said,  there  will  certainly  be  cases  in  the 
field  where  the  non-uniqueness  inherent  to  noisy  inverse  scattering  problems  will  cause  the  whole 
procedure  to  fail  and  yield  dubious  estimates.  In  those  cases  it  will  be  necessary  to  assume  the  target  is 
dangerous  and  dig  it  up. 

In  a  completely  realistic  situation,  where  in  principle  no  training  data  are  given  and  the  ground 
truth  can  be  learned  only  as  the  anomalies  are  excavated,  one  can  never  be  sure  that  the  data  already 
labeled  constitute  a  representative  sample  containing  enough  of  both  hazardous  and  non-hazardous  items. 
This  difficulty  is  mitigated  by  two  facts:  1)  Usually  at  the  outset  we  have  some  idea  of  the  type  of  UXO 
present  in  the  field,  and  2)  The  (usually  great)  majority  of  detected  anomalies  will  not  be  UXO  and  thus 
random  digging  will  produce  a  varied  sampling  of  the  clutter  present.  Methods  involving  semi -supervised 
learning  exploit  this  gradual  revealing  of  the  truth  and  have  been  found  to  perform  better  at  UXO 
discrimination  than  supervised  learning  methods  like  SVM  when  stalling  from  the  point  dipole  model 
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[91,  95].  (Active  learning  methods,  which  try  to  infer  which  anomalies  would  contain  the  most  useful 
information  and  could  thus  serve  to  guide  the  anomaly  unveiling,  show  further,  though  fairly  minor, 
improvement.)  Combining  this  more  powerful  learning  procedure  with  the  excellent  performance  of  the 
HAP/NSMS  method  may  enhance  the  discrimination  protocol  and  should  be  the  subject  of  further 
research. 

In  summary,  the  results  presented  here  show  that  our  search  and  characterization  procedure, 
whose  effectiveness  is  apparent  from  several  recent  studies  [14,  93-94,  96],  can  be  combined  with  an 
SVM  classifier  to  produce  a  UXO  discrimination  system  capable  of  correctly  identifying  dangerous  items 
from  among  munitions-related  debris  and  other  natural  and  artificial  clutter. 

We  repeated  the  analysis  using  the  semi-supervised  Gaussian  mixture  approach.  The  solution 
process  and  results  are  presented  in  Section  6.2.4.  We  found  that  the  method  provides  excellent 
classification  performance  and  has  the  advantage  over  SVM  that  it  is  less  dependent  on  training  data.  This 
made  it  our  preferred  statistical  classification  model,  and  we  have  continued  to  prefer  it. 

6.2.4  Mixed  model  approach  applied  to  Camp  Sibert  data 

We  also  tested  the  mixed  model  approach  of  Section  4.2  on  the  216-sample  Camp  Sibert  data. 
Initially  we  took  the  time  decay  of  the  total  NSMS  over  25  time  channels  for  all  targets  and  parameterized 
it  using  the  Pasion-Oldenburg  law  of  equation  (57).  Taking  the  logarithm  of  that  equation  we  arrive  at  the 
linear  model 

\nQ{t)-\nk-/3\nt-yt.  (166) 

As  features  we  use  k  and  the  ratio  Q(tls)  /  Q(t{) .  Figure  41  is  a  log-log  plot  of  (fir,.  )/  (7(/j  )  vs.  k  . 

Initially  we  used  X-means  clustering  to  estimate  the  number  of  target  types;  the  algorithm  found  five 
clusters  (see  Figure  41).  Then  we  proceeded  to  classify  the  targets.  The  resulting  classification  into  the 
five  classes  is  depicted  in  Figure  42  and  the  corresponding  ROC  curves  are  presented  in  Figure  43. 

The  results  illustrated  that  the  semi-supervised  Gaussian  mixture  model  provides  excellent 
classification  performance  over  the  SVM.  This  made  the  semi-supervised  Gaussian  mixture  our  preferred 
statistical  classification  model,  and  was  used  in  the  consequence  classification  studies. 
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Figure  41:  Log-scale  plot  of  Q(t  )l  Q(tf)  vs.  k  for  Camp  Sibert  data  classification.  Left:  Ground  truth.  Right: 
^f-means  clustering  result. 
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Figure  42:  Classification  of  216  targets  into  five  classes  using  a  bivariate  normal  mixture.  Also  shown  are  the  95% 
confidence  ellipses. 
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Figure  43:  Five  ROC  curves  that  indicate  the  performance  of  the  mixed  model  approach  to  Camp  Sibert  data. 


6.3  Camp  San  Luis  Obispo  (TEMTADS,  MM,  BUD) 

The  discrimination  test  at  Camp  Sibert  UXO  site  was  relatively  simple:  it  involved  discrimination 
of  large  intact  ordnance  from  smaller  clutter  using  data  from  using  first-generation  EMI  sensors.  Real 
sites,  however  contain  assorted  types  of  ordnance,  many  smaller  than  4.2",  and  the  need  to  tackle  this 
more  forbidding  condition  has  prompted  significant  developments  in  both  detection  and  discrimination 
technologies.  Acceptance  of  these  technologies  requires  demonstrating  that  they  can  achieve  100% 
discrimination  confidence  in  terms  of  the  range  of  ordnance  types  and  their  overlap  with  clutter,  while 
taking  into  account  the  terrain/vegetation  at  the  site  and  the  effects  of  the  geological  setting  on  EMI 
sensors  [2,  7,  16,91,94,  97-98], 

To  demonstrate  the  applicability  of  the  classification  technologies  for  a  live -UXO  site  with  more 
challenging  topography  and  a  wider  mix  of  targets -of-interest,  in  2009  ESTCP  conducted  a  second 
discrimination  study  at  the  SLO  live  UXO  site  in  California.  Magnetometers  and  first-generation  EMI 
sensors  were  deployed  to  the  site  and  used  in  survey  mode.  Then  the  BUD  and  TEMTADS  systems  were 
used  to  perform  cued  interrogation  of  the  detected  anomalies.  Simultaneously,  the  MetalMapper  was  used 
in  both  survey  and  cued  modes.  The  collected  data  were  preprocessed  by  data  collection  demonstrators, 
who  performed  background  subtraction,  drift  correction,  and  sensor  positioning. 

The  classification  demonstrators  were  provided  with  calibration  data  sets  for  algorithm  testing 
and  classification  performance  analysis.  The  goal  was  not  only  to  identify  if  the  target  was  harmful,  but 
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also  to  classify  it  completely;  i.e.,  to  identify  its  type,  size,  and  caliber.  The  blind  data  sets  contained  one 
or  more  buried  objects  that  could  be  either  one  of  four  ordnance  items  used  at  the  site — 60-mm  mortar 
shells,  2.76"  rockets,  81-mm  projectiles,  and  4.2"  mortars — or  a  piece  of  clutter.  The  clutter  items  found 
on  the  site  are  UXO  explosion  byproducts  like  partial  mortars  (i.e.,  stretched-out  half-shells),  smaller 
shrapnel,  and  man-made  metallic  clutter;  some  examples  appear  in  Figure  44. 

This  section  presents  the  discrimination  studies  carried  out  on  1282  TEMTADS  and  1407 
MetalMapper  cued  blind  data  sets.  The  total  parameterized  NSMS  amplitudes  were  used  to  discriminate 
TOI  from  metallic  clutter  and  to  classify  the  different  hazardous  objects.  First  we  used  the  combined 
NSMS/DE  algorithm  to  determine  the  total  NSMS  for  each  TOI  from  the  training  data  provided  by 
SERDP.  We  used  the  HAP  method  and  a  combined  dipole/Gauss-Newton  approach  to  validate  the 
location  and  orientation  estimates  given  by  NSMS/DE.  We  then  used  the  inverted  total  NSMS  to  extract 
time-decay  classification  features  for  all  cases  and  input  these  to  several  multi-class  statistical 
classification  procedures  to  perform  discrimination.  Once  our  inversion  and  classification  algorithms  were 
tested  on  calibration  data  we  repeated  the  procedure  on  the  blind  data  sets.  The  inverted  targets  were 
ranked  by  target  ID  and  submitted  to  SERDP  for  independent  scoring. 


Figure  44:  Found  Clutter  Items  on  SLO  UXO  live  sites. 
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6.3.1  The  total  NSMS  for  discrimination 

The  reader  may  recall  from  Section  2.3.6  that  the  initial  amplitude  and  the  decay  rate  of  the  total 
NSMS  depend  on  the  size,  the  geometry,  and  the  material  composition  of  the  object  it  represents.  Early¬ 
time  responses  are  associated  with  surface  eddy  currents  and  the  associated  early-time  NSMS  is  directly 
proportional  to  the  object’s  surface;  at  later  times  the  currents  diffuse  gradually  into  the  object  and  the 
response  is  related  to  the  target’s  volume.  Thus  a  small  and  thin  target  like  the  partial  2.36"  rocket  (Figure 
45)  has  a  relatively  small  initial  NSMS  that  decays  quickly,  while  a  large  object  like  the  4.2"  mortar  of 
Figure  46  has  a  strong  immediate  response  that  decays  slowly,  particularly  along  its  axis  of  symmetry. 

These  considerations  may  be  put  on  a  more  quantitative  footing  through  discrimination  features 
that  summarize  these  characteristics  for  the  different  NSMS  curves.  To  that  end  we  employ  the  Pasion- 
Oldenburg  law  in  its  parameterized  form  (see  Section  2.3.7).  We  fried  different  combinations  of  Baa 

B  ,  and  y  for  discrimination  and  in  the  end  settled  for  log [M  (L.)/  M  (7  )|  and  M  (t,)  as 
features  for  use  with  the  model-based  supervised  clustering  of  Section  4.3.1  (Figure  47). 

6.3.2  SLO  discrimination  results 

The  SERDP  Program  Office  provided  us  with  188  TEMTADS  calibration  data  sets  for  the 
inversion  and  classification  algorithms  testing  performance  analysis.  Our  objective  here  was  not  only  to 
identify  if  a  given  target  was  a  UXO  or  not,  but  also  to  classify  it  completely;  i.e.,  to  identify  TOI  type, 
size,  and  caliber.  We  had  the  same  number  of  calibration  data  sets  for  the  MetalMapper  sensor,  but  we 
used  only  two  data  sets  for  each  TOI,  for  a  total  of  ten  data  sets.  The  blind  data  sets  contained  a  single  or 
multiple  buried  objects  that  could  be  either  one  or  more  TOI. 

We  used  the  188  TEMTADS  calibration  data  to  build  a  catalog  of  expected  total  NSMS  values 
that  were  then  tested  on  the  1282  other  cells.  The  TEMTADS  took  data  over  115  channels  that  span  in 
approximately  logarithmic  fashion  a  lapse  of  time  between  100  ps  and  24  ms.  The  TEMTADS  was 
always  placed  30  cm  above  the  ground.  For  each  data  set  we  run  the  combined  NSMS-DE  and  NSMS- 
HAP  method  [3]  to  determine  object  locations. 

The  target  response  was  approximated  with  set  of  magnetic  dipoles  distributed  on  a  spherical 
surface  of  radius  5  cm.  This  sphere  is  divided  into  17  subsurfaces,  each  of  which  is  assumed  to  contain  a 
magnetic -dipole  distribution  of  constant  density.  Once  the  location  of  the  sphere’s  center  is  determined 
then  the  magnitude  of  each  responding  source  is  obtained  and  the  total  NSMS  is  calculated.  The  inverted 
total  NSMS  curves  for  SLO  TEMTADS  calibration  (green  lines)  and  blind  data  sets  (red  lines)  are 
depicted  in  Figure  45  and  Figure  46  for  partial  2.36"  rockets,  4.2"  mortars,  81 -mm  projectiles,  2.36” 
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rockets,  and  60-mm  mortars.  The  results  indicate  that  the  inverted  and  calibration  total  NSMS  time  decay 
curves  are  similar  and  are  good  discriminators.  Also,  as  the  size  of  the  TOI  decreases  the  inverted  total 
NSMS  time  decay  curves  show  a  larger  spread,  making  them  more  difficult  harder  to  discriminate. 


Time  |m  Sec)  Time  Im  Sec) 

Figure  46:  Inverted  total  NSMS  time  decay  profiles  for  4.2"  mortars  (top  left),  81 -mm  projectiles  (top  right),  2.36" 
rockets  (bottom  left),  and  60-mm  mortars  (bottom  right)  in  the  SLO  TEMTADS  test.  The  green  lines  depict 
calibration  data  and  the  red  lines  correspond  to  blind  data  sets. 


Figure  45:  Inverted  total  NSMS  time  decay  profiles  for  the  2.36"  partial  rocket.  The  green  lines  depict  calibration 
data  and  the  red  lines  correspond  to  blind  SLO  TEMTADS  data  sets. 
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Figure  47:  Result  of  the  supervised  clustering  classification  for  the  SLO-TEMTADS  anomalies  using  the  logarithms 
of  M aa(t j)  and  M  (t  )/M  (t  ).  The  supervised  clustering  has  been  trained  with  calibration  data.  The  red 
markers  correspond  to  clutters  and  the  white  ones  to  TOI. 

We  also  determine  the  Pasion-Oldenburg  parameters  k  ,  B  ,  and  y  for  each  anomaly  from 

equation  (166);  the  inverted  parameters  were  used  in  the  supervised  clustering  algorithm.  We  have 
previously  found  [12,14]  that  the  ratio  of  the  inverted  total  NSMS  at  the  82nd  time  channel  to  that  at  the 
first  time  channel,  which  involves  a  fixed  superposition  of  /?  and  y  ,  shows  discernible  clustering  for  this 
particular  data  set  when  combined  with  the  third  parameter.  The  values  of  log10(  Maa(t x)/  Maa(tg0) )  versus 
logic/  M/r/(t^))  are  plotted  in  Figure  47  for  all  TEMTADS  data  sets.  We  see  that  the  inverted  parameters 

are  well  clustered,  and  for  the  most  part  noticeably  distinct  from  those  of  the  others,  suggesting  that  this 
two-dimensional  feature  space  is  good  for  classification  purposes.  This  suggestion  is  confirmed  by  the 
classification  results  that  appear  in  Figure  48. 
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DartmouthSky  TEMTADS  AdvancedModels  TestSet  TOI 


Figure  48:  ROC  curve  for  SLO  TEMTADS  test  data. 

The  inverted  SLO  TEMTADS  and  MetalMapper  test  data  were  ranked  by  target  type  and  caliber 
and  submitted  to  the  SERDP  office  for  independent  scoring,  which  was  carried  out  by  personnel  from  the 
Institute  for  Defense  Analyses  (IDA).  Our  discrimination  results  are  summarized  in  Figure  49,  Figure  50, 
and  Figure  51.  Our  classification  technique  was  able  to  correctly  identify  all  big  UXO,  (the  2.36",  81 -mm 
and  4.2"  projectiles)  for  both  TEMTADS  and  MetalMapper  data.  The  algorithm  had  only  one  false 
negative  (a  60-mm  mortar)  for  MetalMapper.  In  the  case  of  TEMTADS  the  algorithm  missed  two  2.36" 
rockets  and  five  60-mm  mortars.  These  false  negatives  were  mostly  due  to  small  signal-to-noise  ratios. 
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Figure  49:  ROC  curve  for  SLO  MetalMapper  test  data. 
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Figure  50:  ROC  for  SLO  TEMTADS  data  for  individual  TOI. 
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Figure  51:  ROC  for  SLO  MetalMapper  data  sets:  individual  TOI. 


6.3.3  Comparisons  between  NSMS  and  Dipole  models 

6.3.3. 1  Calibration  SLO-TEMTADS  data 

We  now  compare  the  dipole  and  NSMS  models  as  applied  to  SLO  calibration  data.  The  data  were 
inverted  using  both  gradient  search  and  DE.  For  the  gradient  search  100  initial  guesses  were  used  to  avoid 
local  minima,  with  30  iterations  for  each  initial  guess  to  guarantee  convergence.  For  DE  100  iterations 
were  used.  Once  the  targets’  locations  were  determined  the  dipole  polarizability  matrix  and  the  total 
NSMS  were  determined  and  diagonalized  using  JD.  The  inverted  dipole  tensor  principal  elements  and 
total  NSMS  for  two  calibration  cells  (410  and  489,  shown  in  Figure  52)  appear  in  Figure  53.  The  inverted 
dipole  principal  polarizability  elements  are  seen  to  be  totally  different  for  the  same  60-mm  mortar.  For 
Cell  #489  the  dipole  elements  are  not  symmetric,  and  their  inverted  magnitudes  are  much  higher  than  for 
the  other  cell  even  though  the  targets  and  burial  depths  are  the  same.  The  simple  dipole  model  clearly 
breaks  down  while  the  NSMS  technique  predicts  consistent  results  and  is  more  stable  and  accurate.  It  is 
worth  pointing  out  that  other  researchers  reported  the  same  problem  with  this  cell  when  using  the  dipole 
model  and  overcame  it  using  multiple  dipoles. 
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Figure  52:  60-mm  mortars  actually  found  in  calibration  cells  #410  and  #489. 


Figure  53:  Left:  Principal  elements  of  the  polarizability  tensor  versus  time  for  a  60mm  mortar  in  the  SLO  study. 
Right:  Total  NSMS  time-decay  curves  for  the  same  cases.  The  red  curve  corresponds  to  calibration  Cell  #489  and 
the  blue  curve  to  calibration  Cell  #410. 


6. 3. 3.2  Blind  SLO-TEMTADS  data  sets 

A  similar  performance  was  observed  for  deep  targets  in  blind-test  data.  Figure  54  compares 
library  and  inverted  data  using  the  dipole  and  NSMS  models.  In  this  case  a  60-mm  mortar  was  buried 
35  cm  deep.  Due  to  the  low  signal-to-noise  ratio  the  dipole  model  was  unable  to  predict  stable,  symmetric 
polarizability  tensor  elements,  but  the  total  extracted  NSMS  curves  show  axial  symmetry  and  resemble 
the  60-mm  library  curve  well. 
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Figure  54:  Comparison  between  library  and  inverted  blind  tests  for  the  dipole  model  (left)  and  NSMS  model  (right). 
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Figure  55:  ROC  curves  for  SLO  TEMTADS  and  SLO  MetalMapper  discrimination  studies.  Green  and  red  curves: 
Sky/UBC  dipole  results;  blue  curve:  NSMS  results  obtained  by  our  Dartmouth/Sky  group. 

6. 3. 3. 3  SLO -Discrimination  studies 

Using  NSMS  we  inverted  all  SLO  blind-test  data  sets  and  sorted  them  by  target  ID.  The  same 
anomalies  were  inverted  by  researchers  at  SKY/UBC  using  the  dipole  model.  The  ROC  curves  for  the 
SLO  TEMTADS  and  SLO  MetalMapper  discrimination  studies  are  depicted  in  Figure  55.  The  NSMS 
performs  slightly  better  than  the  dipole  statistical  approach  for  TEMTADS  data.  For  the  SLO 
MetalMapper  data  sets  the  NSMS  shows  higher  false  positives  in  comparisons  with  the  dipole  model,  but 
overall  it  has  only  one  false  negative,  while  the  dipole  model  had  three  false  negatives. 
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Figure  56:  ROC  curves  for  SLO  BUD  discrimination  studies. 

6.3.4  SLO  BUD  data  inversion  and  classification  studies 

The  combined  NSMS-DE  algorithm  was  applied  to  the  SLO  live  site  BUD  data  sets  (539 
anomalies)  and  targets  intrinsic  (total  NSMS)  and  extrinsic  parameters  were  extracted  for  each 
anomalies.  The  discrimination  features  (size  and  shape  information)  were  extracted  from  the  total  NSMS 
time  decay  history  curve  and  anomalies  were  classified  using  the  provided  69  training  data  set.  In 
addition,  the  library  matching  technique,  that  uses  the  entire  time  decay  history  of  the  total  NSMS,  was 
also  used  for  the  classification.  The  inverted  targets  were  ranked  as  TOI  and  non-TOI  items.  The  ROC  for 
the  SLO  BUD  data  sets  is  shown  on  Figure  56.  The  studies  showed  that  only  two  2.36  inch  rockets  were 
misclassified. 
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Figure  57:  SLO  TEMTADS  test  Cell  #16.  Left:  All  25  eigenvalues  vs.  time.  Right:  Four  highest  eigenvalues  vs 
time.  The  target  response  is  weak  and  mixed  with  the  sensor’s  electronic  and  background  noise. 
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Figure  58:  SFO  TEMTADS  test  Cell  #103.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above-threshold  eigenvalues 
vs.  time.  Only  two  eigenvalues  are  above  the  threshold,  indicating  a  low  signal-to-noise  ratio. 


6.3.5  SLO  retrospective  analysis 

During  the  SLO  test  our  algorithms  missed  five  60-mm  mortars  and  two  2.36"  rockets.  The 
missed  anomalies  were  in  Cells  #16,  103,  241,  441,  444,  748,  and  1285.  Figure  57  through  Figure  63 
present  the  results  for  each  of  these  anomalies,  along  with  our  comments. 
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Figure  59:  SLO  TEMTADS  test  Cell  #241.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above-threshold  eigenvalues 
vs.  time.  There  more  than  three  eigenvalues  above  the  threshold,  which  indicates  that  the  cell  contains  more  than  one 
target.  The  curves  decay  fast,  illustrating  that  the  targets  are  small. 


Figure  60:  SLO  TEMTADS  test  Cell  #441.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above-threshold  eigenvalues 
vs.  time.  There  more  than  three  eigenvalues  above  the  threshold,  indicating  that  the  cell  contained  more  than  one 
target.  The  fast-decaying  curves  illustrate  that  the  targets  have  thin  walls  or  are  small. 
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Figure  61:  SLO  TEMTADS  test  Cell  #444.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above-threshold  eigenvalues 
vs.  time.  There  more  than  three  eigenvalues  above  the  threshold,  indicating  that  the  cell  contained  several  targets.  In 
addition,  the  curves  decay  fast,  illustrating  that  the  targets  are  small. 
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Figure  62:  SLO  TEMTADS  test  Cell  #748.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above-threshold  eigenvalues 
vs.  time.  More  than  three  fast-decaying  above-threshold  eigenvalues  indicate  the  presence  of  several  small  targets. 
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Figure  63:  SLO  TEMTADS  test  Cell  #1285.  Left:  All  25  eigenvalues  vs.  time.  Right:  Above-threshold  eigenvalues 
vs.  time.  Again,  the  eigenvalues  indicate  that  there  are  several  small  targets  in  the  cell. 
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6.4  Camp  Butner 

The  former  Camp  Butner  is  a  40,384-acre  site  located  approximately  15  miles  north  of  Durham 
and  straddling  Durham,  Granville,  and  Person  Counties,  all  in  North  Carolina.  The  War  Department 
acquired  the  property  from  private  landowners  in  1942  for  use  as  a  training  and  cantonment  facility 
during  World  War  II.  The  camp  was  primarily  established  for  the  training  of  infantry  divisions  (including 
the  78th,  89th,  and  4th)  and  miscellaneous  artillery  and  engineering  units  [85],  A  large  variety  of 
munitions  have  been  reported  as  used  at  the  former  Camp  Butner,  including  rifle  grenades,  2.36"  rockets, 
37-mm  and  40-mm  rounds,  81-mm  mortars,  and  105-mm,  155-mm,  and  240-mm  projectiles.  Although 
the  historical  records  are  not  definitive,  it  is  thought  that  the  targets  of  interest  at  the  site  of  the  test  are 
mostly  37-mm  and  105-mm  projectiles;  some  of  the  former  have  a  copper  band,  others  do  not.  The  clutter 
items  found  on  the  site  are  for  the  most  part  UXO  explosion  byproducts  like  partial  mortars  (i.e., 
stretched-out  half-shells),  smaller  shrapnel,  and  man-made  metallic  clutter.  An  initial  surface  clearance 
was  carried  out  on  the  site  prior  to  the  collection  of  digital  geophysical  data.  Then  an  EM6 1  survey  was 
conducted  on  two  100'  x  100'  grids  for  site  characterization.  A  surface  clutter  analysis  and  excavation  of 
one  of  these  100'  x  100'  grids  confirmed  the  identities  of  the  targets  of  interest  (TOI),  provided  an 
indication  of  their  depth  distribution,  and  gave  the  demonstrators  some  information  about  the  clutter 
environment  at  the  site. 

At  a  live  site  such  as  this,  the  ratio  of  clutter  to  TOI  is  such  that  only  a  small  number  of  TOI  may 
be  found  in  a  10-acre  area,  far  from  enough  to  determine  any  demonstrator’s  classification  performance 
with  acceptable  confidence  bounds.  To  avoid  this  problem,  the  site  was  seeded  with  enough  TOI  to  ensure 
reasonable  statistics.  Three  types  of  targets — 37-mm  and  105-mm  projectiles  and  M48  fuze  assemblies — 
were  thus  used.  The  survey  data  for  the  study  were  collected  with  a  line  spacing  of  50  cm.  The  detection 
threshold  was  set  to  detect  all  37-mm  projectiles  at  a  depth  of  30  cm  [85],  which  for  the  EM61-MK2 
carted  survey  corresponds  to  a  threshold  of  5.2  mV  in  the  second  time  gate.  Using  this  detection  threshold 
a  first  anomaly  list  was  produced.  This  list  was  used  as  a  stalling  point  for  two  detailed  cued  surveys 
carried  out  using  TEMTADS  and  the  MetalMapper. 

Our  team  processed  both  data  sets  independently  using  our  advanced  EMI  discrimination 
techniques  and  occasionally  requesting  training  data  to  assist  during  the  classification  stage.  The  main 
objective  of  this  section  is  to  demonstrate  the  discrimination  performance  of  the  ON  VMS  model  [99]  in  a 
live  UXO  site  under  realistic  field  conditions;  the  method  is  combined  with  DE  optimization  (the  two-step 
approach  described  in  Section  3.3)  to  determine  the  locations,  orientations,  and  time -dependent  total 
ONVMS  of  the  subsurface  targets.  The  latter  depends  on  the  intrinsic  properties  of  the  object  in  question 
and  can  be  used  for  discrimination.  To  streamline  the  process  we  employed  JD  to  estimate  the  number  of 
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potential  targets  before  inverting.  To  classify  the  targets  in  the  MM  data  sets  we  performed  semi- 
supervised  Gaussian-mixture  model-based  clustering  on  the  total  ONVMS  in  a  process  similar  to  that 
described  in  Section  4.3.1.  We  now  present  the  results  of  our  discrimination  and  classification  strategies 
when  applied  to  the  Camp  Butner  TEMTADS  and  MM  blind  cued  data  sets.  The  SERDP  office  provided 
us  with  2291  cases  interrogated  with  each  system.  We  divided  our  team  into  two  groups:  One  group 
processed  TEMTADS  data  and  the  other  worked  on  the  MM  sets;  each  group  worked  independently 
using  different  classification  strategies.  Each  team  constructed  a  custom  training  list  (amounting  to  less 
than  5%  of  the  entire  blind  data)  and  requested  the  ground  truth  for  those  anomalies  for  use  during  the 
classification  stage. 


6.4.1  TEMTADS  data  discrimination  strategy  and  classification  results  using  supervised 

clustering 

We  processed  all  the  TEMTADS  data  using  the  JD  and  ONVMS  models.  Initially  we  used  JD  to 
estimate  the  data  quality  and  the  number  of  potential  targets.  The  JD  algorithm  constructs  a  multi-static 
response  matrix  using  TEMTADS  data  and  computes  its  eigenvectors  and  eigenvalues,  the  latter  as  a 
function  of  time.  Studies  show  that  these  eigenvalues  are  intrinsic  properties  of  the  targets  and  that  each 
target  has  at  least  three  eigenvalues  above  the  threshold  (noise  level).  For  example,  Figure  65  shows  the 
eigenvalues  extracted  for  a  105-mm  HE  projectile,  a  105-mm  HEAT  round,  an  M-48  fuze,  and  a  37-mm 
UXO.  As  the  number  of  targets  increases  (as  in  Figure  64  and  the  third  row  of  Figure  65),  so  does  the 
number  of  eigenvalues  above  the  noise  level.  We  thus  examined  the  eigenvalues  versus  time  for  each  case 
and  used  them  to  estimate  the  number  of  targets. 
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Figure  64:  TEMTADS  multi-static  response  matrix  eigenvalues  versus  time  for  some  samples  of  requested 
anomalies. 
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Figure  65:  TEMTADS  multi-static  response  matrix  eigenvalues  versus  time  for  a  105-mm  HE  projectile  and  a  105- 
mm  HEAT  round  (top  row),  an  M-48  Fuze  and  a  37-mm  munition  (center  row),  and  two  clutter  scenarios,  one  with 
two  items  (left)  and  another  with  several  (right)  (third  row). 
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In  addition,  based  on  the  eigenvalues’  time -decay  characteristics  we  built  a  custom  training  list. 
For  the  most  part,  the  list  contained  anomalies  that  had  too  many  above-threshold  eigenvalues,  like  the 
samples  depicted  in  Figure  64.  We  requested  two  batches  of  training  data.  The  first  contained  65 
anomalies,  all  of  which  were  clutter;  some  had  six  eigenvalues  above  the  noise  level,  while  others  had 
several  eigenvalues  mixed  with  the  noise.  The  second  batch  consisted  mostly  of  UXO.  Once  we  had  the 
ground  truth  for  all  75  custom  identified  anomalies  we  proceeded  to  invert  all  TEMTADS  data  sets  using 
a  multi-target  ONVMS  algorithm  combined  with  DE.  We  extracted  the  total  ONVMS  for  every  anomaly. 
Armed  with  the  custom  identified  training  list  and  the  inverted  total  ONVMS  for  each  case  we  created  a 
library  for  M-48  fuzes  and  37-mm  projectiles  without  copper  band.  We  did  not  request  training  data  for 
either  of  the  105-mm  UXO  or  for  the  37-mm  projectile  with  copper  band  because  we  already  had 
TEMTADS  test-stand  data  for  these  targets.  The  JD  and  ONVMS  analysis  clearly  showed  the  presence  of 
those  items  at  the  site.  We  implemented  a  library-matching  technique  in  which  we  quantified  the 
mismatch  in  total  ONVMS  between  library  samples  and  blind  items  and  used  it  to  classify  UXO  and  non- 
UXO  items.  The  inverted  total  ONVMS  for  the  anomalies  that  were  classified  as  105-mm  HE  projectiles, 
105-mm  HEAT  rounds,  M-48  fuzes,  and  37-mm  UXO  with  and  without  a  copper  band  are  depicted  in 
Figure  66  and  Figure  67.  All  the  inverted  total  ONVMS  are  seen  to  cluster  well,  and  each  target  has  a 
total  ONVMS  with  features — such  as  its  amplitude  at  the  first  time  channel,  its  decay  rate,  or  the 
separation  between  the  primary  (blue  lines)  and  secondary  (red  and  green  lines)  components  at  different 
time  channels — that  make  it  amenable  to  identification.  (The  most  difficult  differences  to  discern  were 
between  the  M-48  fuzes  of  Figure  66  and  the  37-mm  projectiles  without  copper  band  of  Figure  67).  These 
features  allowed  us  to  classify  targets  as  UXO  or  clutter  and  also  let  us  sort  the  UXO  by  caliber.  With  this 
knowledge  we  created  a  prioritized  dig  list  that  we  cross-validated  using  the  time-decay  curves  of  the  JD 
eigenvalues. 

The  final  prioritized  dig  list  was  submitted  to  the  Institute  for  Defense  Analyses  (IDA)  for 
independent  scoring.  The  scored  results  were  sent  back  in  the  form  of  a  receiver  operating  characteristic 
(ROC)  curve,  which  we  depict  in  Figure  68.  We  can  see  that  a)  of  the  75  targets  that  were  dug  for 
training,  68  targets  were  not  TOI  (shift  along  x-axis)  and  seven  were  (shift  along  y-axis);  b)  for  95%  TOl 
classification  (the  pink  dot  in  Figure  68)  only  seven  extra  (false  positive)  digs  are  needed;  c)  to  classify  all 
TOI  correctly  (the  light  blue  dot)  only  21  extra  (false  positive)  digs  are  needed;  d)  for  increased 
classification  confidence  the  algorithm  requested  an  additional  thirty  digs  after  all  TOI  had  been  identified 
correctly. 
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Figure  66:  Inverted  total  ONVMS  time-decay  profiles 
munition  and  105 -mm  HEAT  round,  and  (bottom)  M-48  I 
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for  four  Camp  Butner  targets:  (top  row)  105-mm  HE 
uze  and  37-mm  projectile  with  copper  band. 
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Figure  67:  Inverted  total  ONSMS  time  decay  profiles  for  a  37-mm  projectile  without  copper  band. 


Sky  Research,  Inc. 


105 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


Butner  Dartmouth  AdvancedModels  None  TEMTADS  Custom  via  TOI 

160 

140 

120 

(/> 

O) 

9  ioo 

g 

°  80 
a> 
n 
E 

I  60 

40 

20 

0 

0  200  400  600  800  1000  1200  1400  1600  1800  2000 

Number  of  Non-TOI  Digs 

Figure  68:  ROC  curve  for  the  Camp  Butner  TEMTADS  test  data. 

6.4.2  MetalMapper  data  discrimination  strategy  and  classification  results  using 
supervised  clustering 

All  Camp  Butner  MM  data  sets  were  processed  using  a  multi-object  ONVMS/DE  code.  The 
combined  procedure  yields  the  total  ONVMS  for  each  anomaly,  which,  like  the  total  NSMS,  is  intrinsic  to 
the  object  it  represents  and  can  therefore  be  used  for  classification  (See  Section  2.3.6).  As  with  the  total 
NSMS,  early-time  ONVMS  responses  are  associated  with  superficial  eddy  currents  and  thus  directly 
proportional  to  the  size  of  the  object’s  surface,  while  late -time  signals  are  due  to  volumetric  currents  and 
thus  proportional  to  the  target’s  entire  volume. 

These  physics-based  features  were  utilized  in  the  supervised  clustering  algorithm.  We  used  the 
ratio  of  the  inverted  total  ONVMS  at  the  30th  time  channel  to  that  at  the  first.  The  values  of 
log  10[M_(f1)/M_(f30)]  vs.  log|(.| M  (t])\  are  plotted  in  Figure  69  (left)  for  all  Camp  Butner  MM  data 

sets.  We  see  that  the  plotted  quantities  exhibit  a  wide  spread  of  values.  To  use  these  features  for  statistical 
classification,  and  for  determining  clusters  and  a  classification  probability  function,  we  started  by 
dividing  the  scatter  plot  of  Figure  69  (left)  into  subsections.  We  then  applied  the  Gaussian  mixture  model 
to  each  subsection  assuming  that  there  were  five  clusters.  From  the  Gaussian  mixture  model  we  extracted 
the  mean  and  standard  deviations  for  each  cluster  and  built  a  global  classification  probability  function, 
depicted  in  Figure  69  (right)  that  depended  on  the  two  feature  parameters.  The  figure  shows  that  there  are 
55  well-separated  clusters.  We  next  created  a  first  custom  training  dig  list  that  contained  55  anomalies, 
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(i.e.,  one  anomaly  for  each  cluster)  and  requested  the  ground  truth.  The  MM  data  for  each  scenario  were 
inverted  using  the  combined  ONVMS-DE  algorithm  as  though  there  were  one,  two,  or  three  targets 
present,  and  the  resulting  total  ONVMS  amplitudes  were  compared.  Whenever  we  spotted  significant 
differences  we  examined  the  curves  visually  (a  sample  case  is  depicted  in  Figure  70)  and,  based  on  this 
examination,  requested  the  ground  truth  for  an  additional  60  datasets.  Once  we  had  the  ground  truth  for  a 
121  custom  training  data  set,  we  classified  all  targets  as  either  TOI  or  non-TOI  items  using  the  probability 
function  of  Figure  69.  The  classification  based  on  the  supervised  clustering  is  plotted  in  Figure  7 1 :  the  red 
circles  correspond  to  TOI,  and  the  green  dots  to  clutter. 
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Figure  69:  Left:  Scatter  plot  for  all  MM  anomalies  based  on  the  extracted  total  ONVMS.  Right:  Probability  function 
for  all  MM  anomalies. 


Figure  70:  Inverted  magnetic  dipole  polarizability  (left)  and  total  ONVMS  (right)  time-decay  profiles  for  MM 
anomaly  #2504.  The  thin  red  lines  show  a  library  sample,  while  the  thick  blue  and  green  lines  show  the  inversion 
results. 
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Figure  71:  Result  of  the  supervised  clustering  classification  for  the  Camp  Butner  MM  anomalies  using  the 
logarithms  of  M  (tt)  and  M  (tt)/  M  (t30) .  The  supervised  clustering  was  trained  with  calibration  data.  The  red 
markers  correspond  to  clutter  and  the  green  ones  to  TOI. 


We  see  that  the  Gaussian  mixture  model  separates  and  clusters  inverted  parameters  well.  The 
clusters  for  the  TOI  are  noticeably  distinct  from  those  of  the  others,  suggesting  that  this  two-dimensional 
feature  space  is  appropriate  for  sound  classification. 

Using  these  results  we  created  a  prioritized  dig  list  for  the  Camp  Butner  MM  anomalies  and  again 
submitted  the  list  to  the  Institute  for  Defense  Analyses  for  scoring.  Our  classification  results  are 
summarized  in  the  ROC  curve  of  Figure  72.  We  see  that  a)  of  the  121  targets  that  were  dug  for  training, 
120  targets  were  not  TOI  (shift  along  x-axis)  and  one  was  (shift  along  y-axis);  b)  for  95%  TOI 
classification  (pink  dot  in  Figure  72)  eight  extra  (false  positive)  digs  are  needed;  c)  to  classify  all  TOI 
correctly  (light  blue  dot)  only  32  additional  digs  are  needed;  d)  for  increased  classification  confidence  the 
algorithm  requested  33  additional  digs  after  all  the  TOI  were  identified  correctly. 

Our  classification  results  for  both  TEMTADS  and  MM  were  scored  independently  by  the 
Institute  for  Defense  Analyses.  The  scores  we  obtained  reveal  that  our  advanced  models  produce  superb 
classification  in  all  cases.  There  were  no  false  negatives,  and  less  than  5%  of  the  anomalies  had  to  be  dug 
to  achieve  100%  correct  classification.  This  is  the  third  time  our  advanced  EMI  and  statistical  models 
have  shown  successful  classification  performance  on  a  realistic  live-site  blind  test. 
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Butner  Dartmouth  AdvancedModels  None  MetalMapper  Custom  v2  TOI 


Difficult  TOIs: 
#884  (37mm) 
#178  (37mm) 
#1728 (Fuze) 
#1002  (37mm) 
#1000 (Fuze) 
#856  (37mm) 
#773  (37mm) 
#720  (37mm) 


Figure  72:  ROC  curve  for  Camp  Butner  MetalMapper  test  data. 


Time  [msec] 


Figure  73:  Left:  Total  ONVMS  time-decay  curves  for  a  105  mm  projectile  in  the  camp  Butner,  NC  study.  Right: 
Principal  elements  of  the  polarizability  tensor  versus  time  for  the  same  case. 

6.4.3  A  Comparison  between  ONVMS  and  Dipole  model 

To  illustrate  the  ONMVS  superior  classification  performances  over  a  simple  dipole  model,  here 
we  analyze  extracted  dipole  polarizabilities  and  total  ONVMS  for  105  mm  projectile.  The  data  were 
collected  at  camp  Beale,  NC  using  the  TEMTADS  sensor.  The  object's  intrinsic  parameters  were  inverted 
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using  the  ONVMS-DE  and  the  simple  dipole  model-DE  algorithms,  with  100  iterations.  The  target  dipole 
polarizability  matrix  and  the  total  NSMS  were  determined  and  diagonalized  using  JD,  and  are  illustrated 
in  Figure  73.  The  results  show  that  magnitudes  of  extracted  dipole  principal  polarizabilities  versus  times 
are  out  of  orders;  namely,  at  early  time  gates  amplitudes  of  the  primary  polarizability  Mzz  are  less  than  the 
secondary  Mxx  and  tertiary  Myy  polarizabilities,  while  the  TONVMS  magnitudes  have  consistent  orders 
for  all  time  channels  and  provides  good  classification  parameters. 

6.4.4  Camp  Butner  retrospective  study  using  semi-supervised  clustering 

In  this  section  we  discuss  an  automated  classification  approach  that  could  provide  a  rigorous 
framework  for  data  processing  that  maximizes  both  the  sensitivity  and  the  specificity  and  minimizes  any 
other  external  inputs  from  the  user.  Such  a  method  would  ideally  infer  the  patterns  from  the  data  itself  and 
determine  the  minimal  possible  training  data  set  required  to  correctly  identify  trends  and  classify  the  data. 
The  ground  truth  could  then  be  requested  from  the  site  and  further  incorporated  into  the  model.  The 
method  would  internally  update  for  the  newly  available  results  and  then  could  either  request  more  ground 
truth  to  further  refine  the  classification  or  provide  a  final  priority-weighted  dig  list. 

Several  clustering  and  statistical  signal  processing  techniques  have  been  applied  previously  to 
UXO  detection  based  on  magnetometry  and  EMI  data  [7-8,  53],  where  classification  was  based  on 
features  extracted  from  direct  magnetic  field  observations.  In  [100],  for  example,  first-generation  (single 
transmitter/single  receiver)  data  collected  by  the  EM61  sensor  was  treated  as  an  image,  with  no  physical 
models  used  for  classification,  and  the  results  were  compared  to  those  using  a  simple  dipole  model. 
Successful  applications  of  Bayesian  data  fusion,  multivariate  Gaussian  representation  of  the  EMI  signals, 
and  semi-supervised  learning  techniques  have  also  been  reported  [95,  100-102],  Here  we  use  a  simple 
combination  of  hierarchical  clustering  and  probabilistic  classification  approaches  to  perform  unsupervised 
(or  semi -unsupervised)  learning  for  UXO  classification  using  NSMS -extracted  Pasion-Oldenburg 
parameters  from  TEMTADS  Camp  Butner  data.  We  first  use  agglomerative  clustering  in  the  feature 
space  to  split  the  entire  data  set  into  a  finite  number  of  clusters  (which  is  an  external  parameter,  and  was 
assumed  to  range  from  1%  to  5%  of  the  number  of  items  in  the  data  set).  Then  we  request  the  ground 
truth  for  the  anomalies  which  lie  closest  to  the  geometrical  center  of  each  cluster  in  the  feature  space. 
Those  clusters  which  then  happen  to  be  centered  around  a  TOl  are  further  labeled  as  potential  UXO 
clusters  and  used  as  a  basis  to  construct  a  Gaussian  Mixture  model  (fitting  either  one  or  more  multivariate 
Gaussian  distributions  across  the  suspicious  clusters).  Any  other  anomaly  can  then  be  assigned  a 
probability  of  being  a  particular  type  of  UXO  based  on  its  position  in  the  feature  space  relative  to  the 
identified  UXO  clusters.  These  probabilities  can  be  used  to  sort  the  anomalies  and  generate  a  prioritized 
dig  list. 
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As  mentioned  previously,  the  multi-target  magnetic  field  inversion  approach  often  provides  better 
target  localization  and  more  precise  NSMS  decay  law,  thus  yielding  more  reliable  features.  In  a  two- 
target  case,  however,  this  means  having  total  of  a  factor  of  three  data  points  in  feature  space  (one  for 
object  1  alone,  second  for  object  2  alone,  and  third  for  objects  1  and  2  considered  acting  as  a  single 
object).  For  every  physical  anomaly,  therefore,  there  is  a  triplet  of  points  in  feature  space  describing  its 
contents.  If  any  of  these  three  points  is  suspected  to  be  a  UXO,  the  whole  “anomaly”  (the  objects  buried 
at  that  specific  location  at  the  UXO  site)  should  be  treated  as  UXO.  It  turns  out  that,  while  it  is  still 
possible  to  perform  the  data  clustering  in  the  feature  space,  there  is  no  straightforward  way  to  identify 
which  training  data  to  request  and  how  to  interpret  it.  Suppose  a  training  data  point  is  requested  from  a 
certain  cluster,  containing  only  clutter  objects  (perhaps  even  having  outlier  values  as  features).  While  the 
data  point  corresponding  to  this  object  indeed  has  features  peculiar  to  clutter,  it  can  happen  that  this  signal 
is  coming  from  a  triplet  containing  a  UXO,  which  will  be  revealed  in  the  ground  truth.  Since  it  is  very 
difficult  to  determine  which  of  the  points  in  the  triplet  belongs  to  which  of  the  physical  items  (object  1, 
object  2,  or  objects  1  and  2  acting  as  a  single  object),  one  would  have  to  mark  all  three  of  the  clusters  as 
potential  UXO,  immediately  leading  to  a  large  rate  of  false  positives. 

To  get  around  this  issue  we  employ  a  two-step  approach.  In  the  first  step,  the  features  extracted 
from  an  inversion  assuming  only  one  target  (single-target  inversion)  are  clustered  and  the  ground  truth  is 
requested  for  the  data  points  closest  to  the  cluster  centroids.  There  are  several  options  for  clustering, 
which  are  taking  different  criteria  into  the  account.  Some  of  the  options  we  found  useful  for  Camp  Butner 
data  classification  are 

1.  Ward  linkage  criteria  with  Euclidean  distances  [103-106].  The  Ward  technique  is  based  on  the 
minimization  of  the  increase  in  the  total  within-cluster  sum  of  squared  distances  between  the 

members  of  a  cluster  and  its  centroid:  E  —  Xf  ,X  „  II  x  .  -  m,  II2 ,  where  K  is  the  total  number  of 

*= 1  xj<ECk  i  k 

clusters  and  m(  is  the  centroid  of  cluster  Ck  [105].  In  agglomerative  hierarchical  clustering, 
therefore,  only  those  clusters  are  merged  which  cause  the  minimal  increase  in  this  distance. 

2.  Weighted  Pair  Group  Method  Average  (WPGMA)  linkage  based  on  Mahalanobis  distances. 
When  any  two  clusters  are  merged,  WPGMA  [104-105]  uses  a  recursive  approach  to  update  the 
distances  between  already  existing  clusters  and  a  newly  formed  one  by  weighting  the  pairwise- 
average-distances  to  original  merged  clusters  with  respect  to  the  number  of  elements  in  them.  The 
Mahalanobis  distance  [107]  provides  a  way  to  measure  the  separation  of  a  point  from  a  particular 
statistical  distribution  described  by  a  given  covariance  matrix  and,  therefore,  takes  into  the 
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account  the  natural  variation  and  spread  of  physically  different  feature  values  in  their  own 
dimensions. 

The  central  elements  belonging  to  each  of  the  clusters  then  are  probed,  and  those  identified  as 
UXO  are  stored.  In  the  second  step,  two-target-inversion  features  are  clustered  with  the  same  algorithm, 
with  the  number  of  clusters  being  3  times  greater  than  that  in  the  first  step.  After  that,  the  confirmed  UXO 
feature  coordinates  for  a  single-target  inversion  can  be  imported  into  the  multi-target  feature  space,  and 
the  clusters  closest  to  these  points  were  marked  as  potential  UXO.  Such  an  approach  combines  the  ease  of 
training  data  interpretation  of  the  single -target-inversion  case  with  the  rigor  and  accuracy  of  the  multi¬ 
target  inversion.  The  multivariate  Gaussian  Mixture  model  can  then  be  constructed  around  the  identified 
UXO  clusters,  and  the  rest  of  the  anomalies  assigned  a  probability  of  being  UXO. 

The  combined  clustering/Gaussian  mixture  approach  therefore  provides  a  natural  way  to  find 
intrinsic  patterns  in  noisy  feature  data  and  yields  a  convenient  probabilistic  measure  of  class  membership 
for  unknown  items.  It  also  reduces  the  amount  of  required  training  data,  improving  both  classification 
sensitivity  and  specificity. 

6.4.4. 1  Results 

In  this  section  we  apply  the  classification  techniques  described  above  to  the  blind  data  set  from 
former  Camp  Butner.  This  data  set  contains  2291  anomalies,  with  no  initial  training  data  available.  The 
suspected  UXO  types  are  37-mm  and  105-mm  projectiles  and  48-mm  fuzes.  The  goal  of  automated 
classification  process  is  to  minimize  the  involvement  of  human  experts  in  the  learning  process.  By 
delegating  the  routine  tasks  such  as  feature  extraction,  clustering  and  labeling  to  the  software,  it  is 
possible  to  extract  only  the  key  structural  information  from  the  complex  data,  leaving  the  less 
cumbersome  but  more  crucial  tasks,  such  as  decision-making  and  quality  control,  to  human  experts.  With 
no  training  data  available,  only  unsupervised  learning  techniques  have  to  be  used  at  the  first  stage  of  the 
process.  Below  we  report  on  the  progress  of  the  blind  test  classification  study  for  Camp  Butner,  which, 
overall,  consisted  of  the  following  steps: 

(a)  The  features  log  k  ,  b ,  and  g  were  extracted  from  EMI  data  sets  of  all  anomalies,  corresponding 
to  1 -object,  2-object  and  3-object  inversions. 

(b)  Initial  clustering  was  performed,  and,  in  order  to  probe  the  feature  space,  the  ground  truth  was 
requested  for  all  targets  whose  features  were  located  closest  to  their  corresponding  cluster 
centroids  (a  total  of  69  targets). 

(c)  Clusters  containing  at  least  one  UXO  were  identified  and  a  smaller  domain  was  selected  within 
the  feature  space  for  further  interrogation. 
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(d)  Second  clustering  was  performed  within  the  selected  domain,  and  those  targets  with  features 
closest  to  the  corresponding  cluster  centroids  were  probed  for  ground  truth  (26  targets).  The 
clusters  with  at  least  one  identified  UXO  were  marked  as  suspicious. 

(e)  All  targets  whose  features  (based  on  2-object  inversion)  fell  inside  any  of  the  suspicious  clusters 
were  used  to  train  a  3 -component  Gaussian  Mixture  model  classifier  and  score  all  of  the  unknown 
targets. 

(f)  All  targets  with  a  score  greater  than  a  specifically  selected  threshold  value  of  log-likelihood  were 
assumed  to  be  UXO,  and  the  ground  truth  was  requested  for  them  (a  total  of  131  targets,  3  of 
which  had  already  been  requested  previously.  Of  these  targets,  118  were  confirmed  UXO). 

(g)  A  new  3-component  GMM  classifier  was  then  trained  based  on  the  features  from  the  3-object 
EMI  inversion.  All  the  items  were  re-scored  to  correct  for  the  changes  and  adapt  for  new 
information.  Another  20  targets  with  consecutively  decreasing  scores  (stalling  form  a  specific 
low  value)  were  then  selected  for  additional  verification. 

At  this  point ,  if  the  verification  yielded  that  all  these  20  targets  are  clutter,  the  algorithm  would 
stop,  and  the  scored  values  would  be  used  to  produce  a  final  dig  list. 

(h)  Four  out  of  the  20  items  requested  happened  to  be  UXO,  and  the  classification  continued.  NOTE: 
the  ground  truth  for  three  of  these  20  items  had  already  been  requested  in  the  previous  steps,  with 
one  being  a  confirmed  UXO. 

(i)  All  confirmed  UXO  were  separated  into  three  groups  (105-mm,  48-mm,  and  37-mm)  without 
further  discriminating  between  the  differences  within  each  group.  Each  of  the  three  groups  was 
used  to  train  a  separate  1 -component  GMM  classifier,  which  was  then  used  to  score  all  of  the 
targets  with  a  separate  score  for  each  of  the  target  types  (based  on  the  features  from  the  precise  3- 
object  EMI  inversion).  The  ground  truth  was  then  selected  individually  for  each  of  the  object 
types,  based  on  a  certain  threshold  score.  This  step  helped  resolve  the  possible  biases  arising  from 
simultaneous  treatment  of  all  targets. 

(j)  A  total  of  36  items  were  requested  from  a  105-mm  scored  data  set,  with  18  being  already  known; 
174  items  were  requested  from  a  37-mm-scored  data  set,  with  118  being  already  known;  and  53 
items  were  requested  from  a  48-mm-scored  data  set,  with  27  of  them  being  already  known. 

(k)  At  this  stage  a  total  of  322  items  were  requested,  162  of  which  were  UXO. 

(l)  Finally,  a  3 -component  GMM  classifier  was  trained  on  the  confirmed  UXO  and  further  used  to 
score  all  of  the  unknown  targets.  A  specific  threshold  was  then  selected  and  the  final  dig  list 
produced.  100%  of  UXO  were  identified  correctly,  with  a  total  of  295  non-TOI  items  (false 
positives).  Total  number  of  anomalies  in  the  data  set  was  2291 . 
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Figure  74:  ROC  curve  for  Camp  Butner  live  site  classification.  100%  UXO  were  identified  correctly,  with  only  295 
false  positive  rate.  The  total  number  of  anomalies  is  2291.  The  blue  dot  corresponds  to  a  threshold  in  the  dig  list, 
when  the  boundary  between  UXO  and  clutter  was  assumed  after  scoring.  The  cyan  dot  specifies  the  actual  position 
of  this  boundary.  In  ideal  circumstances  the  blue  and  cyan  points  will  coincide.  Performing  extra  digs,  however, 
helps  maintain  better  statistics  and  improve  the  results. 
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Figure  75:  Camp  Butner  single-object  inverted  data  clustering.  Left.  Results  of  weighted-linkage  clustering  using 
Mahalanobis  distances  for  single-object  inverted  EMI  features.  Right :  All  four  identified  UXO  (black)  after  a  second 
clustering  within  a  smaller  domain  (logk  <s[2;8|,  h  g[0.05;2|,  g  e[0.05;2]  )  using  Ward  linkage  and  Euclidean 
distances. 


Figure  75  shows  the  results  of  the  first  two  clustering  processes,  corresponding  to  steps  (a)-(d) 
above.  Only  four  UXO  targets  were  identified  at  this  step:  two  37-mm  (with  features  very  close  to  each 
other  in  Figure  75),  one  48-mm  and  one  105-mm. 

Figure  76  illustrates  the  training  data  used  to  create  a  3 -component  GMM  classifier  in  step  (f)  and 
the  resulting  score  distribution  histogram.  The  external  interaction  from  the  expert  in  this  case  consisted 
in  selecting  a  threshold  value  for  scoring  beyond  which  the  ground  truth  would  be  requested.  We  picked 
the  value  of  log(score)  ~  0.5,  which  resulted  in  the  right  peak  in  the  histogram  being  probed,  and  yielded  a 
high  number  of  1 18  confirmed  UXO  out  of  131  probed  items. 

The  newly  acquired  data  was  then  used  to  re-train  the  GM-classifier  (step  g),  using  the  features 
from  a  precise  3-object  EMI  inversion  set  (note  that,  since  the  3-object  inversion  data  set  provides  a  set  of 


Sky  Research,  Inc. 


114 


January  2012 


UXO  discrimination:  Combining  advanced  EMI  models  and  statistical  signal  processing 


MM- 1572  Final  Report 


seven  decay  curves  per  anomaly,  only  the  features  closest  to  already  identified  UXO  centroids  were 
considered  for  GM  training).  The  results  of  the  updated  GM-based  clustering  are  demonstrated  in  Figure 
77.  The  broadening  of  the  histogram  peak  corresponding  to  UXO  is  observed.  Based  on  the  updated 
histogram,  the  ground  truth  from  additional  20  suspicious  items  was  requested  to  statistically  challenge 
the  classifier.  The  region  was  identified  to  be  corresponding  to  log( score)  in  the  region  between  -6  and  - 
5,  based  on  the  external  input  from  the  expert  (visually  observing  the  isosurfaces  and  how  they 
encompass  the  existing  UXO  clusters,  and  considering  the  spread  of  the  histogram  peak  corresponding  to 
the  UXO.  This  step  can  potentially  be  automated  in  the  future  to  increase  process  efficiency).  At  this 
stage,  if  all  of  the  20  items  were  returned  as  clutter,  the  process  would  stop  and  the  scored  items  would  be 
used  to  create  the  final  dig  list.  However,  it  turned  out  that  4  out  of  20  items  were  UXO,  and  therefore  the 
classifier  had  to  be  updated  once  again  to  ensure  that  all  possible  outliers  were  accounted  for. 

In  order  to  resolve  possible  biases  from  simultaneous  treatment  of  different  types  of  targets,  all 
confirmed  UXO  were  separated  into  three  categories  based  on  their  type  (105-mm,  48-mm  and  37-mm, 
without  further  discriminating  between  the  differences  within  each  group),  and  each  group  was  used  to 
train  a  separate  1 -component  GMM  classifier,  which  was  then  used  to  score  all  of  the  targets  with  a 
separate  score  for  each  target  type  (step  (i)).  The  ground  truth  was  then  selected  individually  for  each  of 
the  object  types,  based  on  threshold  score  values  that  were  identified  visually,  using  external  expert  input, 
as  before.  Figure  78,  Figure  79,  and  Figure  80  present  the  results  obtained  with  these  individual  classifiers 
for  37-mm,  48-mm  and  105-mm  target  clusters  respectively. 

The  ground  truth  obtained  as  a  result  of  steps  (a)  through  (k)  constituted  322  requested  anomalies, 
with  160  of  them  being  confirmed  as  UXO.  At  the  final  stage,  a  3-component  GMM  classifier  was  trained 
on  the  confirmed  UXO  from  the  accumulated  ground  truth,  and  further  used  to  score  all  of  the  unknown 
tai'gcts  (Figure  81).  A  specific  threshold  was  then  selected  manually  and  the  final  dig  list  produced.  As  a 
result,  100%  of  the  UXO  were  identified  correctly  (Figure  74),  with  a  total  of  295  non-TOI  items  (false 
positives)  exposed  in  the  process  (total  number  of  anomalies  in  the  data  set  was  2291). 

The  overall  process  yielded  a  successful  ROC  curve.  In  the  future,  the  combined  clustering  and 
GMM  algorithm  should  be  researched  for  further  improvement  to  automatically  find  ways  for  optimal 
data  clustering,  scoring  and  thresholding.  While  external  inputs  from  an  expert  are  valuable  for  guiding 
the  learning  process,  it  is  desirable  to  minimize  human  involvement,  reducing  it  to  a  system-wide  quality 
check  and  classification  control.  For  example,  an  ideal  learning  mechanism  would  first  utilize  all  possible 
information  contained  in  the  data  before  delegating  the  crucial  decision-making  to  human  experts.  Since 
humans  are  still  better  at  some  tasks  involving  pattern  recognition,  matching  or  classification,  such  a 
combined  framework  may  result  in  overall  improved  performance  and  effective  resource  allocation. 
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Figure  76:  Camp  Butner  clusters  used  to  train  the  first  GM  classifier  and  its  results.  Left.  Assumed  UXO  clusters 
used  to  generate  the  3-component  GM  classifier.  Right.  Score  histogram  showing  the  number  of  anomalies  scored 
within  a  particular  range  of  the  loglprobability  density)  in  arbitrary  units.  The  ground  truth  was  requested  based  on 
thresholding  the  log(score)  at  the  externally  selected  value  of  -0.5. 


Figure  77:  Updated  GMM  classifier  after  confirming  118  UXO  in  the  Camp  Butner  data.  Left.  GM-classifier  score 
iso-surfaces  in  the  classification  case  based  on  all  currently  identified  UXO  targets  (118  items),  in  the  feature  space 
corresponding  to  3-object  EMI  inversion.  Right.  Updated  score  histogram  showing  the  number  of  anomalies  scored 
within  a  particular  range  of  the  log(probability  density)  in  arbitrary  units.  An  additional  20  items  were  requested  for 
statistics  to  probe  the  region  corresponding  to  log(score)  within  [-6;  -5]  (this  region  was  identified  with  external 
input  from  an  expert  by  observing  the  corresponding  score  iso-surfaces  and  the  histogram  behavior). 


Figure  78:  GMM  classifier  results  for  37-mm  targets.  Left:  1-component  GM-classifier  score  iso-surfaces  in  the 
classification  case  based  solely  on  identified  37-mm  UXO  targets,  in  the  feature  space  corresponding  to  3-object 
EMI  inversion.  Right :  Score  histogram  showing  the  number  of  anomalies  scored  within  a  particular  range  of  the 
loglprobability  density)  in  arbitrary  units.  A  total  of  174  anomalies  were  requested  (with  118  being  already  known) 
based  on  the  log(score)  cut-off  value  of  about  -6  (specified  externally  by  an  expert  to  allow  enough  statistics). 
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Figure  79:  GMM  classifier  for  48-mm  targets.  Left.  1-component  GM-classifter  score  iso-surfaces  in  the 
classification  case  based  solely  on  identified  48-mm  UXO  targets,  in  the  feature  space  corresponding  to  3-object 
EMI  inversion.  Right :  Score  histogram  showing  the  number  of  anomalies  scored  within  a  particular  range  of  the 
loglprobability  density)  in  arbitrary  units.  A  total  of  53  anomalies  were  requested  (with  27  being  already  known) 
based  on  the  log(score)  cut-off  value  of  about  -5  (specified  externally  by  an  expert  to  allow  enough  statistics). 


Figure  80:  GMM  classifier  for  105-mm  targets.  Left.  1-component  GM-classifier  score  iso-surfaces  in  the 
classification  case  based  solely  on  identified  105-mm  UXO  targets,  in  the  feature  space  corresponding  to  3-object 
EMI  inversion.  Right:  Score  histogram  showing  the  number  of  anomalies  scored  within  a  particular  range  of  the 
loglprobability  density)  in  arbitrary  units.  A  total  of  36  anomalies  were  requested  (with  18  being  already  known) 
based  on  the  log(score)  cut-off  value  was  about  -20  (specified  externally  by  an  expert  to  allow  enough  statistics). 


Figure  81:  Final  GMM  classifier  for  Camp  Butner.  Left:  3-component  GM-classifier  score  iso-surfaces  in  the 
classification  case  based  on  all  identified  UXO  targets,  in  the  feature  space  corresponding  to  3-object  EMI  inversion. 
Right:  Score  histogram  showing  the  number  of  anomalies  scored  within  a  particular  range  of  the  loglprobability 
density)  in  arbitrary  units.  A  total  of  377  anomalies  were  scored  as  UXO  based  on  the  log(score)  cut-off  value  of 
about  -10  (specified  externally  by  an  expert  to  allow  enough  statistics  and  iso-surface  separation  from  identified 
UXO  clusters). 
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6.4. 4.2  remarks 

A  hierarchical  agglomerative  clustering  approach  followed  by  Gaussian  mixture  probabilistic 
modeling  was  applied  in  a  blind -test  format  to  a  live  UXO  site  Camp  Butner.  The  ground  truth  for  a  total 
of  322  items  was  requested  in  a  5-level  iterative  prediction-correction  process,  resulting  in  160  correctly 
identified  UXO.  A  probabilistic  GM  model  was  then  used  for  final  scoring.  Overall,  this  method  yielded 
100%  accuracy  in  UXO  detection,  at  a  cost  of  295-object  false  alarm  rate  (with  a  total  number  of  buried 
anomalies  of  2291).  Machine -learning  techniques  therefore  hold  promise  for  high-quality  automated 
UXO  discrimination,  reducing  the  expert  workload  and  improving  the  process  speed.  Novel  ways  of 
process  improvement  can  be  studied  in  the  future  to  reduce  the  false-alarm  rate  and  improve  overall 
classification  quality.  An  attractive  direction  for  further  research  is  the  creation  of  a  UXO  library 
containing  the  features  extracted  from  UXO  EMI  curves  from  various  live  camps,  and  using  this 
knowledge  in  the  process  of  a  new  live  site  UXO  identification  in  an  automatic  format,  with  minimal 
involvement  by  human  experts. 
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7  Conclusions 

Project  MM- 1572  supported  the  development  of  several  innovative,  robust,  and  noise-tolerant 
EMI  forward  models  and  statistical  signal  processing  methodologies  for  use  in  subsurface  target 
localization,  characterization,  and  classification  at  live -UXO  sites.  In  this  report  we  have  outlined  the 
mathematical  fundamentals,  physical  meaning,  and  practical  realization  of  forward  models  such  as  the 
normalized  surface  magnetic  source  (NSMS)  model  (using  both  charges  and  dipoles  as  sources)  and  the 
orthonormalized  volume  magnetic  source  (ONVMS)  technique.  Both  of  these  procedures  have  been  seen 
to  provide  an  accurate  representation  of  the  EMI  responses  of  subsurface  metallic  targets.  The  models 
were  combined  with  data-inversion  approaches — gradient  search,  direct  search/differential  evolution,  and 
the  like — to  invert  data  collected  by  current  advanced  EMI  sensors.  We  also  developed  and  used  the  HAP 
method  for  estimating  target  locations  directly.  In  addition,  we  explored  several  advanced  statistical 
signal  processing  and  classification  approaches — support  vector  machines,  Gaussian  mixture  models, 
etc. — as  possible  tools  for  discriminating  UXO  from  non-hazardous  anomalies. 

We  adapted  every  model  we  developed  to  a  complete  suite  of  next-generation  sensors,  including 
the  MetalMapper,  TEMTADS,  MPV,  and  BUD.  Comparison  between  gradient  search,  DE,  and  HAP 
showed  DE  to  be  the  most  robust,  noise -tolerant  and  reliable  method  to  determine  extrinsic  parameters  of 
targets;  the  procedure,  moreover,  requires  no  regularization,  and  works  quite  well  when  confronted  with 
multi-target  cases.  For  these  reasons  we  consider  DE  to  be  our  foremost  choice  to  estimate  target  location 
and  orientation.  The  combination  of  DE  with  the  NSMS  and  ONVMS  models  was  extensively  tested  on 
actual  data  and  provided  excellent  agreement  with  the  ground  truth  at  every  instance,  regardless  of  the 
number  of  targets  in  the  cell.  The  models  were  further  combined  with  state-of-the-art  classification 
algorithms  and  applied  to  live-UXO  sites. 

Initially,  we  tested  the  NSMS-HAP-SVM  and  NSMS-HAP-Gaussian  combinations  on  EM-63 
data  taken  by  ESTCP  over  216  test  cells  at  Camp  Sibert  in  Alabama.  The  Gaussian  mixture  model 
provided  excellent  classification  performance,  with  neither  false  positives  nor  false  negatives,  while  SVM 
had  a  tiny  number  of  false  alarms.  In  the  next  test  we  applied  the  NSMS-HAP  and  NSMS-DE 
combinations  to  TEMTADS  data  taken  at  the  APG  standardized  test  site.  We  found  that  the  inverted 
classification  feature  parameters  (the  total  NSMS  in  this  case)  were  well-constrained  for  all  objects  and 
that  the  locations  inverted  using  DE  were  in  good  agreement  with  the  ground  truth.  There  were  214 
anomalies  and  six  types  of  targets  in  the  APG  data  set:  25-mm,  37-mm,  60-mm,  81-mm,  and  (two  kinds 
of)  105-mm  projectiles.  For  each  cell  we  determined  the  total  NSMS,  extracted  discrimination  features 
from  the  NSMS  decay  curves,  and  classified  the  features  using  the  Gaussian  mixture  model  and  a  library  - 
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matching  technique  with  the  help  of  test-stand  and  calibration  data.  The  results  of  independent  scoring 
were  the  following:  1)  All  UXO  were  correctly  identified  as  such  and  correctly  identified  by  type/caliber. 
2)  There  was  a  false  positive  rate  of  ~5%. 

The  classification  abilities  of  the  NSMS-HAP  and  NSMS-DE  algorithms  in  combination  with  the 
Gaussian  mixture  model  and  library  matching  were  again  put  to  the  test  with  data  taken  at  Camp  San-Luis 
Obispo  in  California  using  TEMATDS,  MM,  and  BUD.  There  were  four  types  of  TOI:  60-mm,  81 -mm, 
2.36",  and  4.2"  munitions.  Comparisons  between  the  different  methods  demonstrated  NSMS-DE  to  be 
more  robust  and  stable  than  NSMS-HAP  when  extracting  extrinsic  parameters  from  to  actual  live-site 
data  sets,  particularly  in  multi-target  cases.  This  made  us  adopt  DE  as  our  “official”  procedure  for  target 
pinpointing.  The  blind  test  at  SLO  showed  that  NSMS-DE  can  be  combined  with  the  Gaussian  mixture 
model  and  library  matching  to  reliably  classify  single  well-separated  targets  and  anomalies  with  high 
SNR.  However,  the  method  was  unable  to  identify  all  targets  correctly  (it  missed  respectively  one,  five, 
and  one  targets  for  MM,  TEMATDS,  and  BUD).  We  then  conducted  a  retrospective  study  that  clearly 
demonstrated  the  main  difficulties  at  the  SLO  site:  a  low  SNR  and  the  abundance  of  multi-target  cases.  To 
address  those  issues  we  extended  the  NSMS  technique,  developed  the  ONVMS  model,  and  adapted  the 
JD  method  to  next -generation  sensors. 

The  ONVMS  model  assumes  that  measured  secondary  fields  are  due  to  a  volume  distribution  of 
interacting  magnetic  dipoles;  the  corresponding  Green  functions  are  Gram-Schmidt  orthonormalized  to 
avoid  the  ill-conditioning  and  instabilities  that  plague  multi-object  inversion  and  to  make  the  method  run 
faster.  The  JD  technique,  based  on  diagonalizing  a  multi-static  response  matrix  and  associating  the 
number  of  eigenvalues  above  a  certain  threshold  with  the  number  of  illuminated  targets,  is  reliable  and 
robust  and,  since  it  requires  no  inversion,  essentially  instantaneous.  Additionally,  the  eigenvalues  allow 
one  to  perform  a  preliminary  target  discrimination. 

The  resulting  ONVMS-DE-JD  combined  technique  was  first  used  to  conduct  a  retrospective 
analysis  of  the  SLO  data.  After  that  we  applied  the  procedure  to  yet  another  ESTCP  blind  test,  this  one 
held  at  Camp  Butner,  North  Carolina,  using  the  MetalMapper  and  TEMTADS  instruments.  The 
TEMATDS  and  MM  data  were  analyzed  independently  of  each  other.  The  total  time -dependent  ONVMS 
was  extracted,  inverted,  and  classified  for  each  cell  using  ONVMS-DE-JD  and  both  the  Gaussian  mixture 
model  and  library  matching.  Our  results,  scored  by  the  Institute  for  Defense  Analyses,  consistently 
demonstrated  that  our  methods  do  a  superb  job  of  classifying  anomalies.  There  were  no  false  negatives, 
and  less  than  5%  of  the  anomalies  had  to  be  dug  to  achieve  100%  correct  classification.  A  high-quality 
automated  UXO  discrimination  process  based  on  machine -learning  techniques  has  been  demonstrated  for 
reducing  the  expert  workload  and  improving  the  process  speed. 
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Both  the  SLO  retrospective  study  and  the  Camp  Butner  blind  test  clearly  demonstrated  that  the 
suite  of  advanced  modeling  and  classification  tools  developed  by  our  group  are  robust  and  noise-tolerant 
and  provide  excellent  classification  results  using  real-world  data  collected  by  next-generation  EMI 
sensors.  ONVMS  proved  superior  to  NSMS  and  simple  dipole  model  for  inversion  and  classification 
purposes  and  shall  remain  our  preferred  method  of  analysis.  The  ONVMS -DE-JD  combination, 
supplemented  by  our  classification  algorithms,  was  further  tested  under  ESTCP  Project  201101  using 
MetalMapper,  MPV,  and  2x2  3D  TEMATDS  data  collected  at  Camp  Beale  in  California.  Not  only  were 
the  advanced  EMI  models  able  to  classify  all  “easy  seed  UXO  items”,  they  also  managed  to  identify  all 
other  targets,  no  matter  how  unexpected  or  site-specific,  and  as  small  as  3-cm  fuzes  [108]. 
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